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Abstract 


Several organic compounds were evaluated using fluorescence spectro- 
photometry as possible standards for the determination of trace amounts of 
crude oils or residual fuel oils in sea water. Chrysene was selected as 
the most suitable standard. 


Dichloromethane was evaluated as a solvent for the extraction of oil 
in water at concentrations of a few micrograms per litre and was found 
suitable for the purpose. 


Samples of sea water were collected on Cruise 73-006 of the weather- 
ship, CCGS Quadra. A portion of each sample was analyzed aboard ship with 
the remainder being retained for repeat analysis in the shore-based labora- 
tory in Victoria. The average concentration of fluorescent extractable com- 
pounds in sea water collected at Ocean Station P (50°N, 145°W) was 0.038 ug/1 
when expressed as an equivalent concentration of chrysene. The average 
concentration, expressed as an equivalent concentration of a typical crude 
oil, would be about 0.4 ug/l. 
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Introduction 


A method based on the work of Levy (1971) and Keiser and Gordon (1973) 
was adopted for the determination of low concentrations of petroleum oil in 
sea water. In this method oil is extracted from sea water with an organic 
solvent and its concentration in the sea water is determined by measuring 
the intensity of fluorescence of the extracted oil in a suitable organic 
solvent. The method has several attractive features. It is relatively 
simple to learn and can be carried out quickly, making it suitable for 
routine analysis. It can be used for both quantitative and qualitative 
work with the results being only minimally affected by weathering (Coakly, 
1973). Because of the high sensitivity of fluorescence spectrophotometers, 

1 to 2 litre samples’ of even the cleanest natural waters are sufficient for 
analysis, and because of the high selectivity inherent in the fluorescence 
technique, interferences tend to be less serious than they are in some other 
methods. A particular disadvantage of the fluorescence method, as herein 
described, is that only members of certain small groups of compounds present 
in oils are measurable the bulk of the compounds in oils being undetectable. 
Thus, a small measurable fraction of an oil is being used as a tracer for 
the whole and unfortunately, situations do occur where the fate of the tracer 
is different from the fate of the whole. Consequently, results must be 
interpreted cautiously. 


The measurable groups of compounds in oils have in common the ability 
to fluoresce; that is, to emit light of longer wavelength than the incident 
or exciting light. Fluorescence takes place, moreover, generally within 
microseconds of excitation. Phosphorescence, a related phenomonen, is 
generally characterized by longer time intervals between excitation and 
emission. 


The polycyclic aromatic hydrocarbons (PAH) constitute a major group of 
fluorescent compounds in petroleum oils. Many PAH and some related hetero- 
atom containing polycyclic aromatic compounds are carcinogenic. Teratogenic 
and/or mutagenic properties may also be displayed. Consequently, the fluo- 
rescence measurement of PAH and their hetero-atom containing relatives has 
a value which transcends its value in locating and tracing probable oil 
spills. 


The spectra, emission or excitation, obtained from a fluorescent com- 
pound are characteristic of its electronic structure. Compounds having a 
different electronic structure give different spectra. Since oils differ 
in the relative concentrations of fluorescent compounds, the spectra 
obtained for individual oils, like the spectra obtained for individual 
compounds, will vary, each oil having distinct emission and excitation 
spectra. Since oils also vary in the concentration of their fluorescent 
compounds relative to their non-fluorescent compounds, the fluorescence 
method is inherently most sensitive for the detection of oils which have the 
largest relative concentration of fluorescent compounds. 


The work described in this report was undertaken in part to familiarize 
ourselves with the fluorescent characteristics of different types of oils 
and to assess the applicability of the fluorescence method of analysis of 
some typical oils in sea water. We also sought a readily available organic 
compound as a calibration and reference standard so that other laboratories 
could compare their results with ours. After we had become familiar with 
the fluorescence method and had found a suitable standard, during August 4 
to September 19, 1973, we undertook a study of sea water on Cruise 73-006 
(Fig. 1) of the weathership CCGS Quadra. 


Selection of a Standard 


There are many advantages in using a fluorescent compound rather than 
a fluorescent oil as a reference and calibration standard. A fluorescent 
compound can be chosen which is obtainable from many sources in unlimited 
supply. Purification to meet the criteria for a standard does not repre- 
sent a serious problem. Once it is pure, a fluorescent compound need only 
be protected against chemical degredation, i.e., oxidation, photolysis, 
etc., to guarantee homogeneity between and within batches. In contrast, a 
standard fluorescent oil may be in limited supply from only one source, 
such as the American Petroleum Institute. Individual batches may vary 
unless strict attention to homogeneity is observed. Even within a batch 
homogeneity may be lost if strict handling precautions are not adhered to. 


Since we felt that the oils most likely to be encountered would be 
crude oils, fuel oils or bunker oils, we looked for a standard for these 
oils rather than one for lighter refined oils. Since heavy oils generally 
show excitation maxima in the region 290-330 nm and emission maxima in the 
region 360-400 nm (Appendix 1), these characteristics were taken into con- 
sideration in the selection of the standard compound. Nearly two dozen 
compounds were examined before a choice was made (Appendix 2). For the 
selection of the trial compounds, a table, ''Fluorescence Excitation and 
Fluorescence Emission Spectra of Model Compounds", compiled by McKay and 
Latham (1972), was very helpful and is reproduced in Appendix 3. There were, 
however, some discrepancies between their spectra and our spectra which 
could be attributed in part to the instruments used. Our instrument, a 
Perkin-Elmer Model 204, does not have the resolution of the more expensive 
Perkin-Elmer Model MPF-2A. Variation in detector response shape and source 
characteristics may also have contributed. It should be noted that the 
wavelengths of excitation peaks of solutions containing concentrations of 
over 250 ug/l were not even close to those reported by McKay and Latham. 
Dilution of these solutions resulted in more favorable comparisons. 
Intermolecular solute interactions or interference from undissolved solute 
particles were likely causes of the shift in the excitation peaks of the 
concentrated (possibly saturated) solutions. 


The compound selected for our reference and calibration standard was 
chrysene. Although it is only slightly soluble in cyclohexane and hexane, 
it is otherwise particularly well suited as it has numerous peaks (Fig. 2) 
within the required ranges. Chrysene has excitation maxima at 308 nm and 
318 nm (308 nm and 322 nm according to McKay and Latham) and emission maxima 
at 366 nm, 383 nm, and 400 nm (363 nm, 375 nm, 383 nm and 404 nm according 
to McKay and Latham). Six pairs of wavelengths could be used for comparison 


with heavy oils. The 308 nm, 383 nm pair seemed the best choice since both 
wavelengths lie in the middle of the desired ranges. Further, our assign- 
ment of these wavelengths as the positions of the maxima agrees with that 
of McKay and Latham. The plot (Fig. 3) of emission intensity against 
chrysene concentration (in cyclohexane) shows good linearity over the con- 
centration range, 10-50 ug/1, which corresponds to the concentration range, 
in chrysene equivalents, of the extracts of the oil-spiked water samples 
which we prepared for study. Good linearity was similarly achieved in the 
range 1-10 ug/1, which corresponds to the concentration range of extracts 
of deep ocean sea water samples which are presumably unpolluted. 


Oil and Chrysene Comparison 


Chrysene was compared to available oils ranging from heavy fuel oils 
to outboard motor oils (Table 1 and Table 2). Since wavelengths used were 
more typical of the heavier oils than the lighter oils, the former gave 
much higher emission intensities. For example, crude and bunker oils 
studied exhibited intensities which were 20 to 80 times greater than that 
of Gulf outboard oil and 10 to 70 times greater than that of Esso diesel 
oil. The light fuel oils show excitation maxima between 280 nm and 290 nm 
and emission maxima between 310 nm and 330 nm. Chrysene comparisons, 
therefore, reflect only the tailing effects of the lighter oils. As chry- 
sene was chosen as a standard for heavy oils, another standard would have 
to be found if the lighter oils were to be examined. Where possible 
appropriate spectra of the extracted oil should be used to identify the 
type of oil being measured and determine the credibility of the standard. 


Selection of Solvents 


The decision to use dichloromethane as the extracting solvent was 
influenced by results obtained by Keiser and Gordon (1973) and verified 
in our own laboratory. Because dichloromethane is more dense than sea 
water, multiple extractions using a separatory funnel are easier to carry 
out than they would be with a lighter-than-water solvent such as hexane. 
Because it has a low boiling point (40°C), the danger of losing significant 
amounts of fluorescent compounds during evaporation is minimized. Recover- 
ies of oil from spiked sea water are excellent, being in the range of 
90-100% under most conditions (Table 3, Keiser and Gordon, 1973). One 
drawback, however, is dichloromethane's solubility in water (2% in fresh 
water). Larger volumes of dichloromethane must therefore be used in the 
first extraction to compensate for its solubility. 


Cyclohexane was used to redissolve the residue after evaporation. 
It was chosen over pentane, hexane, benzene and dichloromethane. Cyclo- 
hexane has the highest boiling point, minimizing concentration changes due 
to evaporation. It also yielded the highest emission values for spiked 
samples. Benzene has similar characteristics in both cases, but was 
rejected because itis a strong absorber of short wavelength ultraviolet 
light. Cyclohexane proved satisfactory as a storage medium since dissolved 
oils showed little or no loss of fluorescence even after several months in 
the dark. 


Recovery of Oil from Water 


The recoveries of different oils from waters of various origins were 
determined (Table 3). Recoveries averaging 92%, for example, were obtained 
for three different crude oils which had been added to previously extracted 
surface water from Station P. The loss of fluorescent oil components during 
roto-evaporation of dichloromethane and during redissolution of the residue 
in cyclohexane was determined for two oils (Table 4). Since this loss was 
approximately the same as the overall procedural loss observed for the crude 
oils in Table 3, the recovery of the crude oils in the extraction step of 
our procedure appears to have been nearly quantitative. The recovery of 
refined oil, however, was not quantitative. The reason for the poor 
recovery of the Essomotorboat oil was not pursued in the study. In this 
regard, however, it may be noted that the water used, which was surface 
water obtained from alongside the dock adjacent to the laboratories in 
Victoria's Inner Harbour, may have contained components, such as surfactants, 
which interfered with the extraction efficiency of the dichloromethane. 
Because of our experience with the extraction efficiency of Esso motorboat 
oil from harbour surface water, we were alerted to the necessity of testing 
water samples from a study area for the extractability of the pollutant oil, 
when it is known, in order that concentrations, which are determined by the 
fluorescence method, can be corrected for extraction efficiency. 


Our findings generally agree with those of Keiser and Gordon (1971) 
who reported average percent recoveries for a No. 6 fuel oil of 99% for the 
roto-evaporation step and 90% for the extraction step. We found, however, 
that losses of oil were relatively greater in the roto-evaporation step 
than in the extraction step. This lack of agreement may reflect differences 
in the water and oils used in the two studies. In addition, there were 
procedural differences. We used a higher bath temperature for roto- 
evaporation of the solvent (55-60°C vs. 30°C) which may have accounted for 
our greater loss on roto-evaporation. For the extraction we used a 50 ml 
portion of dichloromethane followed by a 25 ml portion rather than two 40 ml 
portions. It should also be pointed out, the above arguments not with- 
standing, that Keiser and Gordon reported standard deviations of 11% and 7% 
for the roto-evaporation and the extraction steps, respectively. The dif- 
ferences between their results and our results, which come from a small 
number of experiments, may therefore be more apparent than real. 


Field Trial 


The fluorescence method, as it had evolved by August, 1973, was given 
a field trial on weathership Cruise No. 73-006. Sea water samples were 
collected from the port bow with the ship steaming at approximately 5 knots. 
A solvent washed (distilled water, methanol, and dichloromethane, 
respectively) five-gallon stainless steel bucket suspended on a nylon line 
was swung forward at 45° to the ship. At the appropriate instant, determined 
with practice, the line was allowed toslacken momentarily so that the bucket 
dropped mouth first into the water. The first two grabs were used as rinses 
with the third being retained as the sample. (Current practice, however, 
involves taking the first grab as sample. The reason for this procedural 


with heavy oils. The 308 nm, 383 nm pair seemed the best choice since both 
wavelengths lie in the middle of the desired ranges. Further, our assign- 
ment of these wavelengths as the positions of the maxima agrees with that 
of McKay and Latham. The plot (Fig. 3) of emission intensity against 
chrysene concentration (in cyclohexane) shows good linearity over the con- 
centration range, 10-50 ug/1, which corresponds to the concentration range, 
in chrysene equivalents, of the extracts of the oil-spiked water samples 
which we prepared for study. Good linearity was similarly achieved in the 
range 1-10 yg/1, which corresponds to the concentration range of extracts 
of deep ocean sea water samples which are presumably unpolluted. 


Oil and Chrysene Comparison 


Chrysene was compared to available oils ranging from heavy fuel oils 
to outboard motor oils (Table 1 and Table 2). Since wavelengths used were 
more typical of the heavier oils than the lighter oils, the former gave 
much higher emission intensities. For example, crude and bunker oils 
studied exhibited intensities which were 20 to 80 times greater than that 
of Gulf outboard oil and 10 to 70 times greater than that of Esso diesel 
oil. The light fuel oils show excitation maxima between 280 nm and 290 nm 
and emission maxima between 310 nm and 330 nm. Chrysene comparisons, 
therefore, reflect only the tailing effects of the lighter oils. As chry- 
sene was chosen as a standard for heavy oils, another standard would have 
to be found if the lighter oils were to be examined. Where possible 
appropriate spectra of the extracted oil should be used to identify the 
type of oil being measured and determine the credibility of the standard. 


Selection of Solvents 


The decision to use dichloromethane as the extracting solvent was 
influenced by results obtained by Keiser and Gordon (1973) and verified 
in our own laboratory. Because dichloromethane is more dense than sea 
water, multiple extractions using a separatory funnel are easier to carry 
out than they would be with a lighter-than-water solvent such as hexane. 
Because it has a low boiling point (40°C), the danger of losing significant 
amounts of fluorescent compounds during evaporation is minimized. Recover- 
ies of oil from spiked sea water are excellent, being in the range of 
90-100% under most conditions (Table 3, Keiser and Gordon, 1973). One 
drawback, however, is dichloromethane's solubility in water (2% in fresh 
water). Larger volumes of dichloromethane must therefore be used in the 
first extraction to compensate for its solubility. 


Cyclohexane was used to redissolve the residue after evaporation. 
It was chosen over pentane, hexane, benzene and dichloromethane. Cyclo- 
hexane has the highest boiling point, minimizing concentration changes due 
to evaporation. It also yielded the highest emission values for spiked 
samples. Benzene has similar characteristics in both cases, but was 
rejected because itis a strong absorber of short wavelength ultraviolet 
light. Cyclohexane proved satisfactory as a storage medium since dissolved 
oils showed little or no loss of fluorescence even after several months in 
the dark. 


Recovery of Oil from Water 


The recoveries of different oils from waters of various origins were 
determined (Table 3). Recoveries averaging 92%, for example, were obtained 
for three different crude oils which had been added to previously extracted 
surface water from Station P. The loss of fluorescent oil components during 
roto-evaporation of dichloromethane and during redissolution of the residue 
in cyclohexane was determined for two oils (Table 4). Since this loss was 
approximately the same as the overall procedural loss observed for the crude 
oils in Table 3, the recovery of the crude oils in the extraction step of 
our procedure appears to have been nearly quantitative. The recovery of 
refined oil, however, was not quantitative. The reason for the poor 
recovery of the Essomotorboat oil was not pursued in the study. In this 
regard, however, it may be noted that the water used, which was surface 
water obtained from alongside the dock adjacent to the laboratories in 
Victoria's Inner Harbour, may have contained components, such as surfactants, 
which interfered with the extraction efficiency of the dichloromethane. 
Because of our experience with the extraction efficiency of Esso motorboat 
oil from harbour surface water, we were alerted to the necessity of testing 
water samples from a study area for the extractability of the pollutant oil, 
when it is known, in order that concentrations, which are determined by the 
fluorescence method, can be corrected for extraction efficiency. 


Our findings generally agree with those of Keiser and Gordon (1971) 
who reported average percent recoveries for a No. 6 fuel oil of 99% for the 
roto-evaporation step and 90% for the extraction step. We found, however, 
that losses of oil were relatively greater in the roto-evaporation step 
than in the extraction step. This lack of agreement may reflect differences 
in the water and oils used in the two studies. In addition, there were 
procedural differences. We used a higher bath temperature for roto- 
evaporation of the solvent (55-60°C vs. 30°C) which may have accounted for 
our greater loss on roto-evaporation. For the extraction we used a 50 ml 
portion of dichloromethane followed by a 25 ml portion rather than two 40 ml 
portions. It should also be pointed out, the above arguments not with- 
standing, that Keiser and Gordon reported standard deviations of 11% and 7% 
for the roto-evaporation and the extraction steps, respectively. The dif- 
ferences between their results and our results, which come from a small 
number of experiments, may therefore be more apparent than real. 


Field Trial 


The fluorescence method, as it had evolved by August, 1973, was given 
a field trial on weathership Cruise No. 73-006. Sea water samples were 
collected from the port bow with the ship steaming at approximately 5 knots. 
A solvent washed (distilled water, methanol, and dichloromethane, 
respectively) five-gallon stainless steel bucket suspended on a nylon line 
was swung forward at 45° to the ship. At the appropriate instant, determined 
with practice, the line was allowed toslacken momentarily so that the bucket 
dropped mouth first into the water. The first two grabs were used as rinses 
with the third being retained as the sample. (Current practice, however, 
involves taking the first grab as sample. The reason for this procedural 


change is to minimize coating the bucket with material from the ocean's 
surface microlayer.) A subsample of 1.5 litres was transferred to a 

2 litre separatory funnel and extracted with 65 ml of dichloromethane by 
vigorously shaking the funnel for two minutes. After allowing the layers 
to separate for approximately five minutes, the dichloromethane was drained 
into a 250 ml round bottomed flask. A second extraction using 25 ml of 
dichloromethane was then performed and the extract added to the first. 

The combined extracts were then evaporated to dryness using a water bath 

at 55°-60°C aboard ship and a roto-evaporator under reduced pressure with 

a similar bath at 55°-60°C at our shore laboratory. Care was taken to 
remove the evaporating flasks from the water bath immediately after evapora- 
tion of the solvent. The residue, rarely visible, was then redissolved in 
15.0 ml of cyclohexane and compared with chrysene standards using the P.E. 
Model 204 fluorescence spectrophotometer. Appropriate blank corrections 
were made. Reagent blanks were prepared by evaporating 90 ml of dichloro- 
methane from the same batch used for the extractions and taking up the 
residue in 15.0 ml of cyclohexane. All samples, standards, and blanks were 
compared at six different sets of wavelengths using 308 nm and 318 nm for 
excitation and 366 nm, 383 nm, and 400 nm for emission. 


The average value in chrysene equivalents for all samples collected 
and analyzed (Table 5) at Station P was 0.038 ug/1 with a standard deviation 
of 0.017 ug/1. Line P samples (Table 6) excluding Station 5 which was 
obviously contaminated had an average of 0.040 ug/l with a standard 
deviation of 0.025 ug/1. Of the Line P samples, those from Stations 1 and 2 
had comparatively high concentrations. Stations 1 and 2 were over the 
continental shelf and were also stations of more than normal shipboard 
activity. 


Using the sample (conc. 0.034 ug/1) collected on September 12 as being 
representative of the samples taken at Station P, the equivalent concentra- 
tions (Table 7) of five crude oils and five bunker oils for six wavelength 
pairs were calculated. The equivalent concentrations vary from oil to oil 
with the equivalent concentrations of the crude oils being about double 
those of the bunker or residual fuel oils. This result is a reflection of 
the greater relative quantities of fluorescent compounds in bunker or 
residual fuel oils. The amount of variation of the equivalent concentration 
of a given oil from one wavelength pair to another provides a measure of 
the similarity of the composition of the fluorescent components of the oil 
with the composition of the fluorescent components extracted from the sea 
water. The smaller the variation is, the greater the similarity. No varia- 
tion would indicate that the oil in question was present, unweathered or 
otherwise changed,in the sea water. If the relative standard deviation 
(Table 7) is used as a measure of this variation, then, of the oils studied, 
B.C. light crude (0.16 R.S.D.) would be most like the sea water extract. 
With the exception of the Stewart bunker oil, moreover, the sea water 
extract resembles the crude oils more than it does the bunker oils. 
Inspection of Table 7 indicates that the equivalent concentration of a given 
oil may change in a systematic fashion with a change in excitation and 
emission wavelengths. A general decrease in equivalent oil concentration 
appears to be linked with an increase in excitation and emission wavelengths. 
This apparent relationship can be demonstrated more clearly if the equiva- 
lent oil concentration is plotted against a function which is the square 


root of the sum of the squares of the excitation. and emission wavelengths 
for which the value of the equivalent oil concentration was determined. 
This function can be regarded as a radius (Appendix 4) and for the purposes 
of this discussion is designated as r. In Figure 4 is shown the value of r 
for one pair of excitation and emission wavelengths (A ex = 308 nm and 

A em = 366 nm, respectively). Two plots of equivalent oil concentration 
against r are given in Figure 5. The plots show that the equivalent con- 
centration of the "average" crude oil and the "average" bunker oil first 
decreases and then levels off or increases slightly as r increases. 


Although the plots in Figure 5 are informative, a more useful com- 
parison can be made if the equivalent concentration determined at each 
pair of wavelengths for a given oil is first divided by the average of 
all the equivalent concentrations. The ratios which are obtained in this 
manner can be plotted against r as shown in Figure 6. Should the equiva- 
lent oil concentrations determined at each wavelength pair be the same, 
then each of the ratios would be equal to 1 and the ratio against r plot 
would fall on the horizontal straight line shown in Figure 6. The closer 
the oil and sea water extract resemble eath other (at the wavelengths used 
in the measurement of the fluorescence), the closer the curve of the 
equivalent concentration will be to the horizontal straight line. The sea 
water extract can therefore be seen to more closely resemble the "average" 
crude oil than the "average" bunker oil. The standard deviation of the 
ratios from the mean ratio gives a measure of the goodness of fit between 
the curve of ratios and the horizontal line in Figure 6. It can be 
easily shown (Appendix 5) that this standard deviation is in fact the 
relative standard deviation of the equivalent oil concentrations from the 
mean equivalent oil concentration. The relative standard deviations of the 
various oils studied are listed, as noted before, in Table 7. 


In principle, it should be possible to perform the data analysis 
described above for a data set based on a much larger number of pairs of 
wavelengths. For identification purposes an increase in the data base 
would be much prefered. In this case, however, there may be more than 
one ratio determined for a given value of r so that the ratio against r 
plots would give a band of variable width and density rather than a line. 

A simplification can be achieved by using a set of fluorescent measure- 
ments determined for a set of points ( } ex, } em) which lie on a diag- 
onal straight line such as the line shown passing through the origin and 
point (308, 366) in Figure 4. For the above set of points each value of 

r would correspond to only one ratio. It would not be difficult to obtain 
a set of fluorescent measurements for a set of points ( } ex, ) em) 

along the diagonal line as shown in Figure 4. A fluorescence spectropho- 
tometer can be easily modified to scan along such a diagonal line in a 
process called synchronous excitation of fluorescence emission (Lloyd, 197la). 
This process is particularly useful in the analysis of oils and extracts of 
water which contain oil. More useful information regarding oil composi- 
tion can be achieved in a single synchronous scan than can be achieved by 

a single scan with either the excitation or emission wavelength held con- 
stant (Gordon and Keiser, 1974; Lloyd, 1971 a, b, c). 


On August 20, 1973, a series of sea water samples, which were collected 
at three hour intervals over a 24 hour period, were analyzed. The results 


are given in Table 8. The average equivalent concentration for ( A ex, 
XA em) = (308 nm, 383 nm) was 0.058 ug/l with a standard deviation of 
0.021 ug/l. The relative standard deviation was 0.36 which corresponds 
favourably with the relative standard deviation of 0.63 and 0.45 for Line P 
and Station P samples, respectively. If the samplestaken at 1700 hours, 
which contained zooplankton were omitted from the analysis, the average 
would drop to 0.052 ug/1 and the relative standard deviation to 0.26. 
Although further verification would be necessary, the smaller relative 
standard deviation of samples collected hourly compared to that of samples 
collected daily indicates that there may be day-by-day variation in the 
concentration of fluorescent extractable materials at Station P. 


Storage of Samples 


Samples were brought back to the laboratory for analysis in three 
different forms--as sea water, as dichloromethane solutions and as 
cyclohexane solutions (representing three stages in the analytical work-up). 
All samples were stored in glass containers with Teflon lined plastic caps. 
Both the containers and Teflon linings had been previously cleaned with 
dichloromethane. Some caps also had a liner of dichloromethane-rinsed 
aluminum foil covering the interior including the threads. 


Sea Water Samples (Table 9) 


Samples of approximately three litres were preserved with 60 mg of 
mercuric chloride and stored in a freezer aboard ship and a cooler (<5°C) 
at the lab until time of analysis. The extractions performed in the lab 
on these samples yielded slightly lower results; 0.032 + 0.019 ug/1 com- 
pared to 0.038 + 0.017 ug/1 for the shipboard analyses. This difference 
may be due to cleaner experimental conditions in the laboratory rather 
than loss of fluorescene in the stored sample. Even the significance of 
the difference is questionable, however, since it is smaller than the stan- 
dard deviation of either average. 


Dichloromethane 


After extraction of the sea water the dichloromethane was stored in a 
cool, dark place until the time of analysis. Some samples having caps with 
only Teflon liners were definitely contaminated (Table 10) whereas none of 
the samples having caps with added aluminum foil were obviously contaminated. 
The latter samples had an average concentration of 0.046 + 0.016 ug/1 com- 
pared to 0.035 + 0.010 ug/l found for the same sea water samples following 
analysis on shipboard. The results indicate that there was an increase in 
fluorescence on storage. A likely source of contamination would be the 
plastic caps. Simply rinsing the caps with solvent releases greater 
quantities of fluorescent materials than we have encountered in extracts of 
sea water 


Cyclohexane 


Upon completion of the extraction with dichloromethane and subsequent 
evaporation, the residue was taken up in cyclohexane and stored in a dark 
cooler. The consequences of storing extracts in cyclohexane were similar 
to those of storing extracts in dichloromethane. Samples in tubes with 
only Teflon liners became contaminated and those with aluminum foil gave 
mixed results (Table 10). Conclusions on this small amount of data are 
hard to make, but it does seem that meticulous care in preparing storage 
tubes or a more reliable type of sample bottle will be necessary if the 
samples are to be stored in organic solvents. 


Summary 


After an extensive investigation chrysene was chosen as a reference 
and calibration standard for heavy oils, using its excitation maximum at 
308 nm and its emission maximum at 383 nm as the pair of wavelengths to be 
used for comparison. No problems were encountered at sea that would make 
performance of the analyses impractical abcard ship except for the ever 
present possibility of contamination from the ship itself. Storage of 
samples did pose some problems. There may be some loss of fluorescence if 
samples of sea water are stored although our results are not definite 
enough to confirm this. If samples are to be stored in an organic solvent 
medium, a contaminant«free container must be found as many of our samples 
were contaminated. 


Our results indicated that Station P sea water contained an equivalent 
of 0.4 ywg/1 crude oil or 0.2 ug/l bunker oil during the study period. 
This level is considerably lower than levels reported in the Atlantic. 
Levy (1972) published results for concentrations of residual oils observed 
in the northwest Atlantic off the coast of Nova Scotia (January, 1971) 
which averaged between 1.0 and 2.0 yg/1 using bunker C from the ship 
Arrow as his standard. A few months later (April 1971) using the same 
standard, Gordon and Michalik (1971) found the average oil concentration in 
Chedabucto Bay area to be 1.5 ug/l. Further south, Brown, Searl and 
Elliot (1973) found nonvolatile or persistant hydrocarbons to be generally 
between 1 and 12 ug/1 along well travelled tanker routes. Later Keizer and 
Gordon (1973),doing a survey from Halifax to Bermuda,found that 77% of 
their offshore above 200 metre depth samples contained less than 2 yg/1 of 
oil using Venezuelan crude as a calibration standard. It should be noted 
that the samples from the Atlantic were taken closer to shore and in areas 
of higher ship activity than those collected for our study. 


The concentrations we found must only be regarded as estimates of the 
maximum possible quantities of petroleum in sea water since there is no 
evidence to support the supposition that the fluorescent components in the 
extracts are petroleum derived. Even when an extract can be compared to 
standards from a probable oil source, the accuracy may be limited by 
weathering effects or the presence of other fluorescing materials in sea 
water. 
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FIGURES 


Chart showing Line P station positions. 


(a) Scan of fluorescence emission ( \ ex = 308 nm) of a 114 ug/1 
solution of chrysene in cyclohexane and (b) Scan of fluorescence 
excitation ( X em = 383 nm) of 114 ug/l solution of chrysene in 
cyclohexane. 


Linearity of response of Perkin Elmer Model 204 fluorescence spectro- 
photometer. 


Rectangular coordinate system showing the distance r to the point 
(308, 366) and the section in which fluorescent measurements are 
made. 


Plot of equivalent oil concentration against r. 


Plot against r of the ratio of equivalent oil concentration (measured at 
the six pairs of wavelengths) to the average equivalent concentration. 
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Figure 2. (a) Scan of fluorescence emission (A ex = 308 nm) of a 114 ug/l 


solution of chrysene in cyclohexane and (b) Scan of fluorescence 
excitation (A em = 383 nm) of 114 ug/l solution of chrysene in 
cyclohexane. 
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Figure 4. Rectangular coordinate system showing distance r to the point 
(308 nm, 366 nm) and the sector in which fluorescent measure- 
ments are constrained by the requirements of the instrument 
( X ex = 220 nm) and the fluorescence process ( i em = Xex). 
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Figure 5. Plot of equivalent oil concentration against r. (@) 
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TABLES 


The Concentration of Chrysene which has the Same Fluorescent Intensity 


at Selected Pairs of Excitation and Emission Wavelengths as 1 mg/1 of 
a Given Oil 


The Concentration of a Given Oil which has the Same Fluorescent 
Intensity at Selected Pairs of Excitation and Emission Wavelengths 
as 1 mg/1 of Chrysene 


Recovery of Oil from Water using Dichloromethane as Extractant 


Recovery of Oil from Dichloromethane following Roto-evaporation and 
Redissolution in Hexane. 


Equivalent Chrysene Concentrations of Sea Water Samples Analyzed on 
Shipboard at Station P 


Equivalent Chrysene Concentrations of Sea Water Samples Analyzed on 
Shipboard Along Line P 


Equivalent Oil Concentrations of Sea Water Samples Collected on 
September 12, 1973 


Equivalent Chrysene Concentrations of Sea Water Samples During 
24 Hour Time Series on August 20, 1973 


Equivalent Chrysene Concentrations of Sea Water Samples Brought Back 
to Laboratory for Analysis 


‘Comparison of Sample Storage Methods 
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TABLE 1 


The Concentration® of Chrysene which has the Same Fluorescent Intensity at 
Selected Pairs of Excitation and Emission Wavelengths as 1 mg /1°% of a Given 
oil 


ODAex (nO) cf 308 SASF ###308F 318 308fe0, 318 
A em (nm) 366 366 383 383 400 400 
Type of Oil Equivalent Chrysene Conc., ug/l 
Cod liver 070002) 070002. 0.0003." 0.0003 0.0006 -0.0005 
Gulf outboard 0.0087. 0.0066 0.0047 ..0.0038 0.0055 «0.0046 
Esso diesel 0.024 0.015 0.009 0.006 0.007 0.005 
B.C. light crude 0.051 0.045 0.038 0.034 0.060 0.051 
Peace River crude 0.078 0.068 0.063 0.056 0.104 0.091 
Esso crude Oe En) 0.103 077095 0.084 0-163 0.144 
Red Water Gulf crude C5220 0.104 0.100 0.089 0.174 0.149 
Boundary Bay crude O7120 0.109 0.100 0.091 0.164 0-142 
Esso No. 46 0.185 0.168 0.183 0.166 0.366 O2318 
CCGS° Vancouver, bunker 0.185 0.164 0.178 0.163 0.374 0.346 
cccs® Quadra, bunker Q.193 0.173 0.192 0.170 0.376 0. 315 
Esso bunker C 07235 0.205 On217 0.190 0.437 07381 
csst Stewart, bunker 0.279 0.268 0.2351 0,221 0.407 0.389 


@In cyclohexane. 

othe crude oils were samples of pipeline feed oils obtained from the Imperial 
Oil Refinery, Ioco, B.C. The bunker oils obtained from the ships were not 
identified. The other oils were purchased. 


“canadian Coast Guard Ship. 


lees Survey Ship. 
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TABLE 2 


The Concentration of a Given Oil which has the Same Fluorescent Intensity 
at Selected Pairs of Excitation and Emission Wavelengths as 1 ug/l” of 
Chrysene 


A ex (nm) 308 oLe 308 318 308 318 


A_em (nm) 366 366 383 383 400 400 
Type of Oil Equivalent Oil Conc., ug/l 
B.C. light crude 19.6 22.20 20.3 29.4 L6o7 19.6 
Red Water Gulf crude S.3 9.6 4073 ua Dey Ser 6.7 
Boundary Bay crude oo 9.2 10.0 11.0 Oak 720 
Esso crude 8.4 957 LOTS A e9 6.1 6.9 
Peace River crude 126 14.7 Lee Lis6 9.6 Lie 
CCGS Quadra, bunker SY74 So yee 5.9 Zev ae2 
CCGS Vancouver, bunker 5.4 6.1 5.6 6.1 Lat 29 
CSS Stewart, bunker 520 337 4.3 Blas ye 2.6 
Esso bunker C ew) Kaeo Sey 5.2 ri 2.6 
Esso No. 46 5.4 6.0 ys) 6.0 Lisl se | 


“In cyclohexane 
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TABLE 3 


Recovery of Oil from Water using Dichloromethane as Extractant 


Medium - Oil Added Conc., ug/l Ze Recovery” 
Extracted Stn. P Peace River ee Swe: 90, 96 
sea water Esso crude b ey 96, 94 

Red Water Gulf crude 14.0 88 
Victoria Harbour - Boundary Bay maids” 4.0 hele Ie 
sea water, 5m 2 2.0 78 
Boundary Bay crude 4.0 90, 86 
2.0 84 
Distilled water Boundary Bay crude” 2.0 68,. Lu, 
SHPAT SSNS, 
Victoria Harbour Esso outboard© 19.5 55 267 oo 
sea water, surface 39.0 56, 49 


“Recovery after extraction of 1 liter of water with dichloromethane (50 ml, 
25 ml), roto-evaporation and dissolution in “cyclohexane or “hexane. 


Griiwered using Whatman No. 1 filter paper. 


TABLE 4 


Recovery of Oil from Dichloromethane following Roto-evaporation” and Redis- 
solution in Hexane 


O11 Recovery (%) . Average Loss (% 
Esso outboard 0325-90387. 10 
Peace River crude 855. 91 £2 


“Water aspirator, water bath at 55°-60°C. 
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TABLE 5 


Equivalent Chrysene Concentrations of Sea Water Samples Analyzed on Ship- 
board at Station P 


A ex (nm) 308 318 308 318 308 318 

A em (nm) 366 366 383 383 400 400 

3 is Be 
Sample Date Equivalent Chrysene Conc., ug/1 
August 8 0 ei Oo UL 0.080 0.080 0.080 0.11 
August 11 0.075 0.062 0.048 0.041 0.041 0.066 
August 14 0.031 0.024 0.019 0.014 0.014 0.024 
August 17 0.048 0.025 0.031 0.024 0.024 . 0.044 
August 20 0.082 0.067 0.054 0.046 0.046 0.077 
August 23 0.060 0.049 0.040 0.035 0.035 0.062 
August 26 0.076 0.056 0.039 0.031 0.031 0.043 
August 29 0.044 0.037 0.027 0.025 0.025 0.038 
September 1 0.026 OF022 0.016 0.013 0.013 0.026 
September 4 0.037 O03) 0.031 0.025 0.025 0.058 
September 7 0:049°" “0.03915 20-031 10. 025mumr 07025 emma ues 
September 12 06.053 0.044 0.034 0.027 0.027 0.045 


September 15 0.066 0.053 0.045 0.039 0.059 0.074 
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TABLE 6 


Equivalent Chrysene Concentrations of Sea Water Samples Analyzed on Ship- 
board Along Line P 


Seresmcamys) 6). © 2 S06. eieaio. ahr oes0e a raia! =) ROS 2518) | 
A em (nm) 366 366 383 383 400 400 
1973 
Station Sample Date Equivalent Chrysene Conc., yg/1 
2 August 4 node se ned. 0.069 n.d. n.d. n.d. 
) August 4 0.19 0.17 Ogel'7 0.14 Nad . 0.26 
6 August 4 0.048 0.074 OT035 0.031 0.047 0.046 
7 August 4 0.058 0.085 0.038 0.035 0.051 49.051 
ay? August 6 0.045 0.033 0.026 0.022 0.034 0.030 
2 September 16 00032 0.024 Ge 022 0.005 02033 7052 
10 September 17 0.042 0.033 0.028 0.024 0.049 0.046 
9 September 17 0.026 On0ZS 02022 0.019 0.036 0.040 
8 September 17 0.055 0.043 0.034 0.026 {[@05046 0.045 
7 Lee ies 18 0.026 0.020 0.015 0. ORS Q.026 (0.023 
4 September 18 0.076 0.060 0.044 0.039 | Q2060 “0.057 
3 September 18 O°OSL 0.022 0.018 0.015 0.027 0.029 
2 September 19 0.12 0.10 0.089 02085 0.14 Ores 
Al September 19 0.10 0.10 0.085 0.078 O13 Das 
a 


not determined. 
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TABLE 8 


Equivalent Chrysene Concentrations of Sea Water Samples during 24 Hour Time 
Series on August 20, 1973 


SSS a eR FRR PSS PES Sa TIN BR RRSP OE SEIT FD RS SN TR TE EE I! 


A ex (nm) 308 318 308 318 308 318 

} em (nm) 366 366 383 383 400 400 
Cos | Equivalent Chrysene Conc., ug/l 
0200 0.085 0.076 0.060 0.058 0.10 0.086 
0500 0.10 0.081 0.067 0.054 0.10 0.078 
0800 0.095 0.080 0.063 0.050 0.086 0.072 
1100° 0.058 0.045 0.037 0.031 0.052 0.047 
1400 0.082 0.067 0.056 0.046 0.077 0.071 
1700°2° O16 0.12 0.10 0.085 0.12 0.12 
2000” 0.086 0.071 0.051 0.037 0.072 0.052 
23002 0.053 0.040 0.031 0.028 0.041 0.040 


(Paint chips in sample. 


Pstored overnight in test tubes (cyclohexane medium) with Teflon lined 


plastic caps. 


¢Zooplankton in subsample. 
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TABLE 9 


Equivalent Chrysene Concentrations of Sea Water Samples Brought Back to the 
Laboratory for Analysis 


» ex (nm) 308 ic 308 318 308 318 


rs . 

Station Sample Date. Equivalent Chrysene Conc. pg/l 
2 August 4 0.86 0.068 0.069 02057 use me) O.LE 
E August 8 0.12 0.090 CFO 71 0.054 0.10 OF 75 
zg August 11 0.065 0.049 0.035 0.041 0.035 
P August 14 0.050 0.038 C2052 0.026 0.043 0.041 
in August 17 0.051 0.036 0.032 0.022 0.035 (0.0385 
P August 20 0.048 Os032 02027 OSO 22 0.036, @:05L 
p? August 23 0.086 0.072 0.060 07052 On 2d 0.10 
i August 26 0.061 0.046 02035 0.029 0.046 0.040 
Ne August 29 0.033 0.023 0.018 02013 0,021 O.022 
P September 1 O2033 02025 0.017 0.014 GO. O22 0.022 
2 September 1 0.037 O.QE7 0.013 O2010 AOL -O<OL7 

Duplicate 
P September 4 0.030 O0L7 0.013 0.010 0.018 0.014 
e September 4 0.025 0.014 0.011 0.007 0.013 -. 0.610 
Duplicate 
P September 7 0.026 0.019 On OnZ 0.010 0.016 0.027 
S September 12 0.044 0.029 0.025 0.019 0.034 0.032 


Se See ee es. 


“aluminum foil liner corroded - flakes in sample. 
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Appendices 


Fluorescence Excitation and Fluorescence Emission Spectra of Specific 
Oils 


Fluorescence Excitation and Fluorescence Emission Spectra of Model 
Compounds--This Study 


Fluorescence Excitation and Fluorescence Emission Spectra of Model 
Compounds--McKay and Latham 


Transformation from Rectangular to Polar Coordinates 


Relationship between Relative Standard Deviation of Equivalent 
Concentrations and Standard Deviation of Ratios 
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Appendix 3 


a ene UEEdEEEtEIIEESIEESSEE ESE 


\ 
Table I. Kluorescence Excitation and Fluorescence Emission Spectra of Model Compounds" 


Compound 


Pluorene 

Naphthalene 
9-Methylcarbazole 
(N-Methylearbazole) 
Carbazole 
2-Methylearbazole 
3-Methylearbazole 

11 H-Benzo[4]fluorene 
(2,3-Benzofluorene) 

11 7/-Benzo[a]fluorene 
(1,2-Benzofluorene) 
Triphenylene 

11 H-Benzo[a]carbazole 
(1,2-Benzocarbazole) 
7H-Benzo[c]carbazole 
(3.4-Benzocarbazole) 
Chrysene 

Phenanthrene 
7H-Dibenzo[c,g]carbazole 
(3,4,5,6Dibenzocarbazole) 
Dibenz[a,cJanthracenc 
(1,2,3,4-Dibenzanthracene) 
Picene 
(1,2,7,8-Dibenzphenanthrene) 
Anthracene 
13//-Dibenzo[a,/Jcarbazole 
(1,2,7,8-Dibenzocarbazole) 
Pyrene 

Benz[aJanthracene- 
(1,2-Benzanthracene) 
Benzo[b]chrysene 
(2,3,7,8-Dibenzphenanthrene) 
Dibenzo[y.p]chrysene 
(1,2,7,8-Dibenzochrysene) 
Diben2{a,/anthracene 
(1,2,5,6-Dibenzanthracene) 
Naphtho[1,2,3,4-de/Jehrysene 
(1,2,4,5-Dibenzopyrene) 
Benz{cJacephenanthrylene 
(3,4-Benzolluoroanthene) 
Benzo[a]pyrene 
(1,2-Benzopyrene) 
Dibenzof[c,g]phenanthrene 
(3,4,5,6-Dibenzophenanthrene) 
Benzo[y/iJperylene 
(1,12-Benzoperylene) 
Dibenzo[de/,mmo]chrysene 
(Anthanthrene) 
Benzo[rs/pentaphene 
(3,4,9, 10-Dibenzopyrene) 
Perylene 

Coronene 
Dibenzolh,de/Jchrysene 
(3,4,8,9-Dibenzopyrene) 
Fluoranthene 

Ovalene 


Fluorescence excitation 
wavelength, nm 


29 Seas 


288 
322 


268 275(s) 
269 278 
249 293 


249 
250 
252 
270 


293 
297 
296 
288 


320 
319 
325 
306 


319 326 342 


255) 265) 296) 3063 318 


262 
255 


277 
279 


288 
306 
263 286 324 
261 
261 
278 


Dil 
278 
303 


297 308 
285 296 
350 363 


322 


269(s) 279 289 


287 304 328 
260 


290 


SID Waeby SSN) aot Sey) 
3 


5 
321 334 348 


308 
255 


B22 


271 290 317 329 


25S 289 


280 340 353 


301 351 


276 330) 342358 


255 


256 300 33399349 


239 B32 


291 348 385 


260 384 407 


247 2977, B32 
410 
340 


399 422 


2555) 
DoD 
272 


370 
303 
300 


388 
324 
Se 448 
241 
314 


280 
399 


290 311 
422 448 


256 
328 


266 
342 


spectra, 


344 360 


376 


365) 385 


422 430 


Mos) Bes 


326 344 


395 


300 


303 
J24 
334 
335 
335 
342 
342 
347 


354 
354 


362 
363 
348 
367 
HY 
B77, 


380 
380) 


374 
387 


394 
395 
396 
397 
403 
405 
405 
399 
432 
434 
440) 


411 
45] 


409(s) 


4S0) 


« The most intense peak in each spectrum is in italic type. Shoulders are indicated by (s). 


McKay and Latham (1972) 


Fluorescence emission spectra, 


310 
338 
349 
351 
350 
358 
351 
365 


364 
372 


369 
375 
357 
386 
388 
398 


40| 
40) 


379 


450 


466 
422 
480 


418(s) 
462 475 482 


wavelength, nm 


350 
365(s) 378(s) 
364(s) 

365(s) 385(s) 
372(s) 378(s) 
359 369 


382(s) 


373 381(s) 391(s) 
392 413 

381 402 
383 


365 
405 


404 427 
385 


398 409s) 


421 449 


424 
425 


451 


384 
435 


389 395 
462 


444 470(s) 
418 446 473 
420 446 
429 459 489(s) 
429 457 

446(s) 

420 446 

494 

462 480 494 

SO) S40(s) 

ae 435, 446 455° 475° 485 S08 


436 463 


490 503 509 514 539 


Le ee ee See oa ra i ee 
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Appendix 4 
Transformation from Rectangular to Polar Coordinates 


A rectangular coordinate system is shown in Figure 4 with emission 
wavelengths ( A em) marked along the ordinate and with excitation 
wavelengths ( A ex) marked along the abscissa. Each point in the plane 
has rectangular coordinates ( A ex, } em) and represents a single pair 
of wavelengths. The distance r from the origin (0,0) to the point 
( A ex, A em) is given by the Theorem of Pythagoras: 


r 24 A ex2 + i em2 


The distance r may also be regarded as the radius of a circle with centre 
(0,0) passing through the point ( A ex, A em). Whether r is regarded 

as a distance or a radius, however, a given value of r does not correspond 
to a unique wavelength pair or point ( dex, Ad em). For example, in 
Figure 4 the point (308 nm, 366 nm) has a value of r equal to 478 n, but 
this value of r could correspond as well to any other point of the same 
distance from the origin. The point (308 nm, 366 nm) is uniquely determined 
by the value 478 nm and the value and angle -@ shown in Figure 4. This 

value of @-is 49.9° and is determined using the formula: 


em 


-@-= tan ! z 
A ex 


Thus, the point (308 nm, 366 nm) can be also represented as (478 nm, 49.9°) 
which are called the polar coordinates of the point. 


Sy) 


Appendix 5 


Relationship between Relative Standard Deviation of Equivalent Concentrations 
and Standard Deviation of Ratios 


| n 
s = a r (Gaex)s by definition 


atte eh) 


n-1 1. 
x2 
a ca tenet 
li Bes a 
i= xX 
ee i 2 
“yn-1 Me a 2 
i=l 
n 
but R ape SUMERY: 
n ‘ i 
i=l 
n Xx 
= oe = 
i=l x 
n 
mi n eh xi 
x 


x | |n 
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where s = standard deviation of equivalent concentrations 
about the mean equivalent concentration 
n = number of equivalent concentrations 
th 
Xs =i equivalent concentration 
1 2 
— == Y x, = mean equivalent concentration 
x n Zz ' 
i=l 
* = relative standard deviation of the equivalent 
x concentrations about the mean equivalent 
concentration 
x 
R, = se i yatio 
x 


R = mean ratio 


sl= standard deviation of ratios about the mean 
ratio. 


For additional copies or further information please write to: 


Department of Fisheries and the Environment 
Institute of Ocean Sciences, Patricia Bay 
512 - 1230 Government Street 
Victoria, B.C. 
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THE INTERACTION GF CHLORINE AND SEAWATER 


by 


R.W. Macdonald and C.S. Wong 


Institute of Ocean Sciences, Patricia Bay 
Victoria, Buc. 


February, 1977 


This is a manuscript which has received only limited circulation. 
On citing this report in a bibliography, the title should be followed 
by the words ''UNPUBLISHED MANUSCRIPT" which is in accordance with 


accepted bibliographic custom. 


Abstract 


An examination of the interaction of chlorine and seawater has been made 
with particular emphasis on those aspects most pertinent to a spill or leak 
of liquid chlorine in seawater. This has been done in order to be able to 
predict the results of a rupture of the chlorine tank cars lost in 
Malaspina Strait early in 1975. 


The appropriate phase diagrams have been constructed from information 
available in the literature, and small scale release experiments from the 
submersible, Pisces, confirm that the behavior of chlorine when released 
at various depths conforms rather closely to thermodynamic prediction. 
Depending on depth or temperature, four phases are possible: liquid 
chlorine, gaseous chlorine, liquid water (with chlorine dissolved in it), 
and chlorine hydrate (solid). 


Literature data has been used to predict the effect of chlorine on such 
seawater parameters as alkalinity, pH, and Eh, and the speciation which 
should occur with various chlorine in seawater mixtures. 


Results of experiments performed here show that, at low concentrations 
(5 mg chlorine/1 seawater) seawater has a natural demand of about 1.5 mg/1 
in half an hour, and 3 mg/1 after a day or so. Sediments can also destroy 
chlorine with 5 g (wet weight) removing 5 mg of chlorine from a liter flask 
of stirred seawater within four hours. 


Consideration of critical data on chlorine, and a review of our release 
experiments show that gaseous explosion will not occur when liquid chlorine 
and seawater are mixed under conditions possible in Georgia Straits. 


In the event of a chlorine leakage from a tank car resting on the bottom, 
liquid chlorine and chlorine hydrate will tend to fall being denser than 
seawater, and stay near the bottom. At depths less than 40 meters it is 
possible that gaseous chlorine could form. Depending on depth , much of 
the gaseous chlorine would not reach the atmosphere but would be tied up as 
an hydrate and remain in the seawater. Eventually after sufficient dilution, 
the natural chlorine demand of seawater and sediments will remove the 
chlorine through oxidation reactions. Mixing will help this process. 
Secondary toxic products formed by chlorination (or bromination) of organics 
should not be a problem in a clean area such as Malaspina Strait. 
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Introduction 


On February 19, 1975, four tank cars filled with liquid chlorine were 
lost from a barge being towed in the vicinity of Malaspina Strait in B.C. 
coastal waters (Fig. 1). Although a series of searches was carried out 
(A.D. O'Connor, March 1975, November, 1976 and others), the exact location of 
the sunken cars has not as yet been determined. Each tank car carried approx- 
imately 94 metric tons! of liquid chlorine with a 10% void space padded to 
6.9 bars with air. 


Although there has been no evidence to date which might indicate that 
chlorine has leaked from the cars, it is only a matter of time before cor- 
rosion leads to failure of the tanks. How long this will take is an open 
question since it is not known exactly what condition the cars were in when 
they were lost or what sort of damage they might have received as a result of 
the accident. 


The tank cars are constructed of five circumferential sections of 2 cm 
steel with two layers of 5 cm cork for insulation and a0.3cm steel jacket on 
the outside. The steel should not crack on impact, but rather deform. Strain 
introduced in this manner should not significantly accelerate corrosion, and 
the outside of the tank car should degrade rather slowly, of the order of 
.5 mm per year. Corrosion will not take place on the inside of the tank car 
provided water does not leak in, since a protective chloride layer is formed 
by interaction of the chlorine and iron. The cars are apparently quite rugged, 
and it is reported that when two cars fell down an incline east of Prince 
George, one car punctured the other car, but the second car survived the fall 
and was retrieved intact (anonymous). 


The most vulnerable area for corrosion is around the hot rolled monel 
valve seat. Small scratches in a good coat of paint would lead to pit cor- 
rosion in the iron in close contact with the seat. How quickly this could 
cause penetration depends very strongly on how well the valve body and steel 
were coated with paint. A minimum time of two to three years has been 
estimated, while 10-20 years is more likely. Once a pinhole is formed it 
should expand rapidly since the acidic chlorine water mixture is quite cor- 
rosive (Carson, 1975). Since rupture is inevitable it was decided to examine 
the chemical aspects of the chlorine-seawater system in order to predict the 
likely sequence of events following tank car failure. 


The study was initiated by a literature search out of which the appro- 
priate phase diagrams have been constructed. Further, the likely speciation 
of chlorine and its effect on pH and pE have been estimated. In order to 
fill in gaps in our knowledge of the chlorine seawater system, we ran some 
chlorine demand experiments in seawater obtained from Saanich Inlet and on 
surface sediment from the Strait of Georgia. 


Ist units are used throughout this report. Conversion factors to other 
commonly used units are to be found in Table lL. 


Before relying on thermodynamic predictions we thought it would be 
worthwhile to carry out some small scale chlorine releases at various 
depths of the water column to get some visual verification of the behaviour 
of chlorine in seawater. This was done from the submersible, Pisces IV, 
in Saanich inlet: 


Thermodynamics 


The phase diagram for the chlorine-water mixture forms a three dimen- 
sional system with the variables composition, pressure and temperature. 
Figures 2-5 present some two dimensional aspects of the phase diagram which 
cover conditions likely to be encountered in Malaspina Strait. 


In order to calculate the position of the lines on the phase diagrams, 
data was extracted from Ketelaar (1967) for the chlorine rich side, and 
from Fernandez et al (1967) for the remainder of the diagram. Vapour 
pressures of chlorine were drawn from the Matheson Gas Data Book (1966). 
In making the various equilibrium calculations it was assumed that only a 
two-component system existed, that is, seawater (35°/..) and chlorine. 
Therefore, it cannot be rigorously applied to the case where other atmos- 
pheric components such as nitrogen and oxygen are included. If an atmos- 
phere of fixed composition were also included, there would be four dimensions 
on the phase diagram and the possibility of a gas phase even at high 
pressures. This does not limit the utility of these phase diagrams since, 
as a first approximation, chlorine and seawater mixtures below the surface 
are essentially isolated from the atmosphere. With an unlimited volume of 
atmosphere in intimate contact with the chlorine-water system it is expected 
that chlorine would partition into the atmosphere according to its vapour 
pressure. This may be considered as a possible removal path of chlorine 
from seawater. Chemical reactions in the water column will also lead to 
chlorine removal. 


A brief description of how the phase diagrams were calculated is 
included below to make clear the estimates and assumptions which have been 
used. 


In considering water solubility in liquid chlorine, chlorine solubility 
in water and the composition of the hydrate it is assumed that over a range 
of 0-20 bars, hydrostatic pressure does not cause any appreciable effects. 
It is also assumed that hydrostatic pressure does not significantly alter 
vapour pressure. It can be shown that this is justifiable since 


dP 


a. = liquid (1) 


<| I< 


vapour 


where P is the vapour pressure, P_ is the total pressure (including hydro- 
static and other contributions) and the V's are the respective molal volumes. 
For water vapour, increasing the total pressure from 1-10 bars results in 
less than a 1% increase in vapour pressure. 


Vapour pressures for water over seawater (35°/,.) were taken from 
Sverdrup, Johnson and Fleming (1942). The mole fraction of water in the 
saturated gas phase was estimated as 
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Vapour pressure of water over hydrate was estimated from the data of 
Ketelaar (1967). 


Equilibrium between chlorine gas and water was calculated from the 
data of Bozzo in Fernandez et al (1967) on chlorine solubility in pure 
water NaCl solutions at various pressures and temperatures. (A salting-out 
coefficient for 35°/.. seawater based on the NaCl solutions was calculated 
to be about 0.4) Chlorine solubilities in the hydrate region were obtained 
from Bozzo's Figure 7 in Fernandez et al (1967). The solubilities of 
chlorine in water outside the hydrate region at 9 bars pressure were extra- 
polated for various temperatures from Bozzo's Table IV taking into considera- 
tion the salting-out effect. Figure 6 summarizes the solubility of chlorine 
in seawater as a function of pressure and temperature. 


The solubilities of water in liquid chlorine were taken from the data 
provided ty Ketelaar (1967) for both the liquid water and hydrate regions. 
Here it was assumed that the influence of salt on the activity of water was 
negligible. (This assumption is quite acceptablesince the activity of the 
water is directly related to its vapour pressure, and only a relatively 
small vapour pressure lowering is observed on going from distilled water to 
seawater) 


Fernandez et al (1967), using the Miller-Strong method, determined that 
the formula for chlorine hydrate was Cl, > 6.2 H,O and this composition 
has been assumed for the phase diagrams enrcnented flere. Roozeboom (1918) 
gave the upper invariant point for the hydrate (quadruple point) as 28.7°C, 
and a pressure of 6.08 bars. Recent measurements of Fernandez et al (1967) 
indicate that the quadruple point should be at 28.3°C with a pressure of 
about 8.1 bars. The recent data is probably more accurate and so has been 
adopted here. Salt depresses the temperature of the upper invariant point, 
and this effect can be taken into consideration by reference to Bozzo's 
Figure 6 in Fernandez et al (1967). The temperature at which the hydrate 
becomes unstable at 1 bar was also determined from this source. 


Implications of the Phase Diagram 


The phase diagram can be used to predict a possible sequence of events 
if chlorine is injected at some point into the water column. Since the 
phase diagrams are based on equilibrium calculations, conclusions can only 
be drawn as to what is possible, or thermodynamically acceptable. Other 
aspects such as very slow kinetic behaviour may prevent thermodynamic 
equilibrium from being obtained. With this in mind, let us consider two 


possibilities, a spill at the surface (1.01 bar) and a spill at 80 meters 
(9 bars). 


Reference to Figure 2 shows that at 1 bar, pure chlorine can form a 
gas phase over the whol- range of temperatures from 0-30°C. On mixing 
chlorine with water in a mole ratio of 99:1 there will be a possibility 
of obtaining chlorine gas over a chlorine saturated water solution, or 
chlorine gas over a solid hydrate (crystalline yellow) if the temperature 
is below about 8°C. Further dilution with water will not result in any 
phase changes until the molar ratio of water to chlorine exceeds 6.2:1 at 
which point a water phase in equilibrium with the hydrate is possible 
provided the temperature is sufficiently low. Eventually, when the mole 
fraction of water exceeds .998, only a water phase will be present with 
some dissolved chlorine forming a seawater "bleach". This will occur when 
there are 10 g or less of chlorine dissolved in 1 litre of seawater. 


Reference to Figure 3 shows that at 9 bars pressure one does not 
expect a gaseous phase over the 0-30° temperature range. Following the 
dilution path of chlorine, a liquid chlorine phase is anticipated initially. 
With sufficient water the hydrate will start to form (below 26°C) or chlorine 
and liquid water will coexist above that temperature. Further dilution will 
eventually result in dissolution of the hydrate or liquid chlorine. 


If the bulk of seawater at Malaspina Strait is at or near 8°C (Crean 
and Ages, 1968) reference to Figure 4 will give a rough idea of the various 
possibilities. Only in a very small region is a gas phase possible, and 
it is anticipated that initially liquid chlorine and hydrate would occur 
(below 40 m) and after further dilution only hydrate and water would be 
thermodynamically stable. Eventually all hydrate would dissolve given a 
sufficient quantity of water. 


Speciation 


Chlorine reacts with water very quickly to form hypochlorous acid 
which can then further dissociate 


Cl, + H,0 * HOCI + HORE CB (3) 


HOC] 2 HY + ocl™ (4) 


The equilibrium constant for Equation 3 has been measured by at least 
two groups and their results have been reviewed by White (1972). 
AE 255C, 
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The dissociation of hypochlorous acid in seawater has recently been 
examined by Sugam and Helz (1976) and the apparent ionization constant 
Ki has been given as a function of salinity and temperature. 


xt -LOCLEN] 4, (6) 
{ HOC1] 
Where ocl, is the total concentration of hypochlorite, and 4, is the 


operational hydrogen ion activity obtained from a pH measurement with a 
glass electrode and a reference electrode with liquid junction. Figure 7 
shows the fraction ocl, as a function of pH and temperature. 


HOC1 

On a recent cruise on the CSS Vector (November 29 - December 3, 1976) 
a pH profile was determined at two stations. Measurements were made with 
a glass electrode and a reference electrode with liquid junction. The method 
of measurement was based on the experimental design given by Zirino (1975) 
for shipboard determinations. The pH's measured at 25°C have been corrected 
to in situ values by using Table VI.3 in Strickland and Parsons (1972). 
Figure 8 shows the in situ pH as a function of depth, and the station 
locations are marked on Figure 1. As can be seen, the pH of the deeper 
water is at about 7.7, while the pH of the surface water (above 30 m) is 
slightly higher at about 7.9. In seasons where biological productivity is 
higher, some CO, will be removed from the water, and it is anticipated that 
the pH could rise to about 8.3, at least in the surface waters. The tem- 
perature of the deep water varies between 8° and 9°C, while surface waters 
drop to as low as 5°C in winter and rise to 15°C in summer (Crean and Ages, 
1968). The inset area in Figure 7 shows the region normally expected in 
unpolluted water. For sufficiently dilute chlorine the OCl. form will 
dominate with the ocl, ratio varying from 3 to 10. The higher value is to 


HOC1 
to be expected in surface waters particularly when the pH is higher. This is 
of some importance since evidence shows that the non-ionized form is more 
toxic to bacteria (Fair et al, 1968) and also more reactive in the formation 
of toxic secondary products (Morris, 1967 and Farkas et al, 1949). 


There is, however, another factor which should be considered. The 
addition of chlorine to seawater will cause a pH change according to 
Equation 3 and Equation 4 since each mole of chlorine produces initially 
one mole of HCl, and further ionization of hypochlorous acid can potentially 
liberate a second proton. 


Consideration simultaneously of both Equation 3 and 4 allows the 
equilibrium concentration of the chlorine species to be calculated as a 
function of pH. The ratio of each species to the total chlorine present is 
shown in Figure 9 where [c1,] 


(7) 


PO ys [Cl] + [HClO] + [C10, ] 


[HC10 ] 


[c1,] + [HC10] + [C10,] (8) 


oT = 


[c10, ] 


oe [Cl,] + [HClO] + [C10] ae 


Consideration has been taken of the chloride concentration in seawater, 
which tends to shift Equation 3 to the left on going from fresh water to 
seawater. From Figure 9 it is important to note that C10, comprises 1024 
or more of the chlorine species in seawater only when the pH is greater 
than 6. From pH 4 to 6 the HC10 form is most important, while chlorine 
becomes significant below pH 4.5. 


Oxidation reactions involving HC10 or C10 can also liberate protons 
which no longer have an associated chlorine species as shown in Equation 10. 


HOC1 = sHCl 37s 0 (10) 


2 
(There is of course a chloride generated by this reaction). When some of 
the chlorine has been lost by this route, it will appear as if a strong 
acid, HCl, has been added along with the remaining chlorine. 


The Effect of Chlorine on Alkalinity and pH 


According to Equation 3, the addition of chlorine to seawater can be 
considered as a titration of seawater with a strong acid. Hydrogen ions 
can also be added through the secondary ionization of HOC1. (This secon- 
dary reaction will not affect the titration alkalinity since for each mole 
of protons added, one mole of the conjugate base OCl is also added. 
Titration to a pH 4 endpoint will cause all of the OCl to be converted 
back to HOC1. As mentioned above, oxidation reactions removing the OC1l- 
form (Equation 10) will lower the alkalinity) . There may also be a further 
titrational effect dependant on the pH of the end point since some chlorine 


will be formed at lower pH's. Reference to Figure 9 will illustrate this 
clearly. 


In the long term, once all of the active chlorine species have been 
reduced to chloride, two moles of H will have been generated for each mole 
of Cl, added to the seawater. Consider the tank cars which contain a total 
of 376 metric tons of chlorine. If we accept an alkalinity change of 2% 


as being ee eg insignificant, simple calculations show that dilution 
with 2 x 10®m3 of seawater are required. 


Let us now consider the influence of adding Cl, to seawater on the pH 
of seawater. A "typical" sample of seawater might fiave a total CO. of 
about 2.2 x 10 ? mol/1 and a pH of 7.9. Using the apparent dissocfation 
constants for the CO, system (Horne, 1968) and Kw taken from Culberson 
and Pytkowictz (19735 for 8°C and 35°/,. salinity, it is possible to 


calculate what additions of strong acid are required to cause the pH to 
change to any specified value in a closed system (i.e., CO, is not allowed 
to escape to an atmosphere). Figure 10 shows the titration curve calculated 
for C, = 2.2 x 10 3? mol/l. Keeping the initial pH at 7.9 and allowing 
effect shown in Figure 10. 


Titration with chlorine has some similarities to titration with a 
strong acid, but cognizance must be taken of both Equation 3 and 4. Using 
the appropriate equilibrium constants of Sugam and Helz (1976), and White 
(1972), the titration curve has been calculated and is shown on Figure 9 
along with the speciation curves. The interesting features are that at 
high pH (6.5-7.9) the system is buffered to a certain extent by OC1l_ 
particularly at about pH 7.4 when [O0Cl ] = [HOC1]. The pH, therefore, 
does not change as rapidly in this region with the addition of chlorine 
as it does when most of the OC1 is converted to HOC1 below pH 6. Below 
pH 3.5 Cl, becomes an important species, since the equilibrium of Equation 3 
is shifted to the left. Addition of large amounts of chlorine here do not 
have a great influence on pH. This is particularly dramatic since the 
logarithm of chlorine added is plotted as the ordinate. 


The measured pH values in Malaspina Straits are fairly uniform with 
depth (Fig. 8) indicating that at least at these two stations at the time 
of sampling no anomalous behaviour atributable to chlorine was determined. 


Effect of Chlorine on Redox Potential 


The redox potential of a system at equilibrium is expressed as a 
dimensionless quantity pE. Somewhat analogous to pH, pE is the negative 
logarithm of electron activity. It is referenced to a conventional zero 
assigned to the standard hydrogen electrode. 


Eh Eh 


z = eo 
223 ERT 59 et 


pE = -logfe} = 


where F is the Faraday, R the gas constant (8.31 joules/°Cmole and T the 
absolute temperature. The redox potential is normally discussed in terms 
of both pE and Eh, the interrelation being demonstrated in Equation ll. 


For oxic seawater the measured pE is found to be approximately 8.5 
suggesting that pE is controlled by the 0,/H,0, couple (Breck, 1975). 


The addition of chlorine to seawater would be expected to upset this 
precariously poised system, chlorine being a strong oxidant. If ityis 
assumed that the seawater system becomes controlled by the C1,/C1 couple 
then 


1, Ci eo Cl =~=s—«E® = 1.3595 (Horne, 1968) (12) 


From fundamental thermodynamics 
AG° = -nFE° = -RT 2.3 log K (13) 


where AG° is the standard free energy change and K is the equilibrium 
constant associated with the respective reaction. 


For Equation 12, AG° = -13.12 KJ/mol 


Tog loge. 25402 (14) 
log K = log [ soe len = log {Cl } + pE -'s log P (15) 
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For seawater of 35°/.. salinity 


{ere yi bs of er *] ie eo oer cn= merce (16) 
( - was taken from Whitfield, 1976) 
Gi 
23.02 = -.044 + pE - %log Py (17) 
Z 
pE = 23.46 + %log P (18) 


cl, 


Variation of Pol from 0.1 to 10 bars will not cause a large change 
2 

in pE and therefore pE = 23.5 should be sufficiently accurate when the 
couple in Equation 12 controls the redox potential. This indicates in 
more precise terms what we already intuitively suspected; that chlorine 
is strongly oxidizing. Reversing the procedure, it obtains directly 
from Equation 18 that normal oxic seawater with a pE of 8.5 has a 
negligibly small Pol since log Pol Rear 3 8 i! 
2 Z 

Perhaps a more useful way of looking at the redox potential of 
chlorinated seawater is to consider the hypochlorite ion. 


clo + HOW iae ye 20H + Cl (19) 


obs 


0.89 (Horne, 1968) 


AG® = -171.5 = -RT 2.3 log K (20) 


log K = 20.14 = 2 log{OH } + log{Cl } + 2pE - log{H,0} - log{clo } (21) 


S = “ 

= 2[log Ks + pH] + log{Cl } + 2pE - log{H,0} - log{c10 } (22) 
= 26.40 + 2pH - 0.44 + 2pE + .008 - log{c10 } (23) 
(Kew has been taken from Culberson and Pytkowicz (1973), for a temperature 


of 25°C and salinity of 35°/... The activity of water in seawater has been 
estimated from the vapour pressure lowering of 35°/.. seawater, Sverdrup 
et al, 1942). 


Therefore 
pE = 28.49 - pH + slog{c10 } (24) 


Using the approximation given by Sugam and Helz (1976), for the 
activity coefficient of the hypochlorite ion in seawater 


3 + 
OCl = 0.64 = PKGL (25) 


pE = 28.39 - pH + glog[C10, ] (26) 


Equation 26 is simply a different formulation of Equation 18, but it 
is easier to see the influence of chlorine on pE because by reference to 
the titration curve in Figure 9 it is now possible to estimate pH. 


Use of Equation 5 and 6 along with the estimated value for pH yields 
pE as a function of chlorine added. This curve is shown in the top part of 
Figure 9 where it can be seen that pE does not vary greatly over a wide 
range of total chlorine concentrations. This curve was generated by 
assuming that no chlorine was removed through oxidation reactions. 

As indicated earlier, oxidation processes will remove chlorine, but will 
not remove the HCl generated and therefore a variety of time-concentration 
paths are possible along the pH and pE curves. These paths will be more or 
less independent depending on how much oxidizable material there is in 
seawater. 


Equation 18 and 26 are also useful in determining the oxidation state 
of the elements or compounds in seawater at equilibrium, and for instance, 
it can be shown that bromide should be converted to bromine at equilibrium 
in the presence of chlorine. 


Kinetics 


There is general agreement that the adsorption of chlorine by water 
(Equation 3) is extremely fast (Shilov and Solodushenkov, 1936, 
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Morris and Carrell, 1946, Spalding, 1962). In fact, Shilov and Solodushenkov 
(1936) found the reaction to be almost complete in less than one second at 
1°C, a temperature somewhat lower than any expected in Georgia Strait 

waters. Spalding (1962) considered two possible reactions 


4 i " 
Cl. + H,0 < HOC] + H + Gl (3) 


clit OH + HOC] + Cl (27) 


His results indicated that between pH 3 and 10.5 the first equation 
applied, and considering a copious quantity of water is available, Equation 3 
can be kinetically treated as a first order equation. Spalding measured a 
first order rate constant of 20.9/sec at 25°C with an activation energy of 
62.8 KJ/mol. The half life of the reaction at 8°C can be calculated at 
about .1 sec, and therefore conversion of cl, to HOC] is very rapid. 


The production of chlorine hydrate (CL, * 602°H230) °appéars ‘to*be trans- 
port controlled with the slowest step occurring at the chlorine-water inter- 
face (de Graauw and Rutten, 1970). Therefore, agitation or turbulence will 
encourage the production of the hydrate. 


The Chlorine Demand of Seawater 


So far we have considered the various equilibrium reactions of chlorine 
in water assuming that no removal mechanism were operating. This, of course, 
is not the case in nature and a variety of removal paths are available. 

Given sufficient time and mixing, all of the chlorine with a +l valance will 
eventually be removed through redox reactions. A variety of compounds 
naturally present in seawater will be important in the removal of chlorine, 
some of which are listed below (White, 1972). 


Ammonia 

Amino Acids 
Proteins 
Organic Carbon 
Nitrite 

Iron 
Manganese 
Sulphide 


Onan Lf WD NE 


Many of the listed compounds are produced and regulated by biological 
processes, and are contained in dead or living material. The most important 
group of compounds with regard to short term chlorine interaction are those 
containing nitrogen (1, 2, 3 and 5 of the above list), especially ammonia. 
The reaction of chlorine and ammonia in dilute solutions leads to the forma- 
tion of the chloramines. 


aL 


HOC] + NH, + NH,C1 + H,0 (28) 
monochloramine 

HOC1 + NH,Cl + HHC1, + H,0 (29) 
dichloramine 

HOC] + NHC1, > NC1, + HO (30) 


nitrogen trichloride 
(trichloramine) 


The above reactions are not redox reactions in that chlorine still 
maintains its +l valance and these compounds still have oxidizing, disin- 
fecting and biocidal capacity. An excellent and very complete review of 
the interaction of chlorine with ammonia and also other compounds is given 
by White (1972), to which the reader is referred for more detail than is 
provided here. Jolley (1973) has also collected together much data on the 
speciation and exchange reactions involving chlorine in water. 


Chlorine which is tied up by ammonia or nitrogen compounds is known as 
"combined available chlorine residual" or just "combined residual". Residual 
chlorine existing as hypochlorite ions or as hypochlorous acid is termed 
"free available chlorine residual" or "free residual" (ASTM Standards, 1976). 
In addition to the interaction of chlorine with nitrogen compounds, other 
reactions are popoenee with the substances listed above. Reducing agents 
such as H,S, Mn2 » Fe* and organic carbon can be oxidized and so remove 
chlorine quickly from the water. Loss of chlorine in this manner is referred 
to as the "chlorine demand" of the water (Singley, 1973). Chlorine demand, 
pea is that amount of chlorine which the water can remove by reducing the. 
Cite L ~l or irretrievably tying up the chlorine in some organic compound. 
Operationally, demand is defined as the "difference between the amount of 
chlorine applied to a treated supply, and the amount of free combined or 
total available chlorine remaining at the end of a contact period" (APHA, 1971). 


With respect to chlorination, seawater has one important difference from 
fresh water. It contains an appreciable quantity of bromide, with 67.3 mg/Kg 
in seawater of 35°/,. salinity (Wilson, 1975). As mentioned while discussing 
the effect of chlorine on the pE of seawater, it is thermodynamically possible 
for chlorine to displace the bromide in seawater with a loss of about 30 mg/Kg 
of chlorine, and a production of about 67.3 mg/Kg bromine in 35°/o. seawater. 
Bromine is then free to react in an analogous manner to chlorine producing 
hypobromous acid and the bromamines. Most instrumental methods do not 
differentiate bromine and chlorine contributions, and therefore seawater 
measurements have commonly been referred to in terms of the chlorine demand. 

A method which reportedly could differentiate between chloramines and broma- 
mines was devised by Johannesson (1958), but Duursma and Parsi (1976) could 
not obtain satisfactory results in their laboratory. In the remainder of 
this report we shall refer to all measurements in terms of chlorine, tacitly 
realizing that some bromine or brominated species will be contributing. 


lhe 


Griffin (1944) in experimenting with higher than normal residuals used 
for water treatment, discovered the so called "breakpoint" phenomenon. 
Measured chlorine residuals initially increase with the amount of chlorine 
added, but suddenly drop off. This behaviour is caused by the interaction 
of chlorine and ammonia. After the minimun, or "breakpoint" free residual 
chlorine can be measured. 


A breakpoint curve was determined for some 2 m seawater obtained in 
Howe Sound from the ship M.V. Pandora II during an Ocean Chemistry cruise 
in July, 1976. Various amounts of chlorine were added to the seawater and 
the total residual was determined amperometrically after 30 minutes. 
(The amperometric determination of free and combined residual chlorine is 
well documented by the ASTM Handbook (1976) and the Fischer and Porter 
Instruction bulletin for the amperometric titrator made by that company. 
The reader is referred to these for details of the method. The Fischer 
and Porter Model 1771010 Amperometric Titrator was used for all chlorine 
determinations, and does not distinguish between chlorine and bromine contri- 
butions.) Figure 11 shows the essential features of this breakpoint curve 
and several aspects should especially be considered. The curve does not pass 
through the origin due to some initial chlorine demand caused by some species 
in seawater which are readily oxidized within the half hour chlorination 
period. The first part of the curve lies in the region where most of the 
chlorine is tied up as monochloramine. After the maximum, dichloramine may 
be formed, an important compound since it participates in further reactions 
which lead to the eventual oxidation of ammonia by deamination. The downward 
slope in Figure 11 results in a measured increase in chlorine demand. The 
minimum is the breakpoint beyond which free chlorine residual exists. In 
fresh water there is also a small portion of dichloramine and trichloramine 
present. Duursma and Parsi (1976) in some very recent investigations 
measured the interaction of ammonia and chlorine in seawater with relatively 
high doses of both chlorine and ammonia (up to 710 mg/Kg Cl,). Their experi- 
ments point out some fundamental differences between the formation of chlora- 
mines in seawater and in fresh water. It would appear that the breakpoint in 
seawater is obtained at lower ratios of chlorine to ammonia than in fresh 
water, and that the dichloramine does not persist in detectable quantities. 
Oxidation of ammonia via the chloramines can result in various possible nitro- 
gen species including Ny» NO., N,0, and NO (White, 1972). 


While investigating the chlorine demand of water as a function of time, 
Taras (1953, 1950) found that chlorine demand could be expressed by the 
equation 


(31) 


where D is chlorine demand, t is contact time, and k and n are constants 
determined experimentally. 


The exponent n relates to the speed of the chlorine consuming reactions, 
and as such helps to define which species are responsible for the chlorine 
demand observed. The higher the exponent, the slower the reaction, and also 
the more complicated the reacting species. 
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Some experiments were made on the chlorine demand of seawater as a 
function of time on samples obtained in August 1976 from Saanich Inlet in the 
vicinity of the small CEPEX enclosures. Both a 2 meter and a 30 meter 
sample are shown in Figure 12 (each seawater sample was chlorinated with 
5 mg/l chlorine initially). The first chlorine demand measurement appears to 
be low possibly due to experimental error or perhaps due to an induction 
period peculiar to the seawater system. (The 2m sample also had the first 
point somewhat lower than might be expected.) The exponent of t was found to 
be .22 for the 30 m sample and .20 for the 2 m sample, indicating that 
compounds other than ammonia are also responsible for some of the chlorine 
demand. Water with the nitrogen component mainly due to ammonia have 
exponents in the .01-.08 range (White, 1972). A Taras plot (Fig. 13) was 
also made for a sample of surface sediment obtained in Georgia Strait several 
miles off Point Grey. In this case 5.2 grams of sediment (wet weight) were 
placed in 5 litres of demand-free water which was then chlorinated to 5 mg/1. 
The sample was continuously stirred by a magnetic stirrer with a Teflon 
stirring bar throughout the measurements. Figure 12 shows that again, the 
early measurement was depressed, while the exponent of t calculated from the 
later points was .28 indicating the likely occurrence of more complicated 
nitrogenous compounds than ammonia. The 5 g sample also had a fairly high 
demand removing all of the chlorine in the 5 litre flask within four hours. 
(All demand experiments were run under ambient laboratory light conditions. ) 


These experiments show that seawater may be expected to have a short 
term demand of about 1.5 mg/1 (Fig. 11). Sediment appears to have a fairly 
high demand with 1 g (wet weight) of surficial mud capable of quickly removing 
2 mg of chlorine from a 1 litre volume. This agrees well with the estimate 
given by White (1972) for the 10 minute demand of seawater of 1.5 mg/l. 
In any event the chlorine demand should not exceed 3 mg/1 unless some organic 
pollution has taken place. Examination of Figure 12 shows that after two days 
only 3 mg/1 chlorine had been removed from the 30 m seawater sample. The same 
was true for the 2 m sample. 


Duursma and Parsi (1976) investigated the persistance of chlorine in 
seawater using far higher concentrations of chlorine than used in our studies, 
with chlorine doses of over 1600 mg/1. They found that under these circum- 
stances a more or less rapid loss of chlorine was often seen with the higher 
dose during the first few minutes. Aeration augmented this effect causing 
a seawater solution initially dosed with 1400 mg/1 of chlorine to fall to 
about 1100 mg/1 after ten hours. They hypothesized that some of the chlorine 
may have been removed by gas bubbling,an unlikely mechanism since at the 
pH of their chlorinated seawater (greater than 10) the amount of free chlorine 
is very low (see Fig. 9). There does not at present appear to be any reason- 
able explanation for this behaviour, and their investigations at lower doses 
of chlorine (below 600 mg/1) do not show this loss. 


A possible direct removal of chlorine from seawater is possible by a 
light stimulated reaction 


ay 
Cl, + H,0 + 2HC1 + 0, C32) 
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Duursma and Parsi (1976) studied the effect of light on persistence of 
chlorine in seawater with chlorine concentrations below 35 mg/1. The found 
no significant difference in chlorine behaviour between the light and dark 
conditions, and suggest that formation of bromine may be responsible since 
the bromine reaction analogous to Equation 32 was not photochemically 
mediated under the experimental conditions. 


Gaseous Explosion 


When a liquified gas is spilled onto or injected into a second fluid, 
there is a possibility of forming superheated liquid drops. In this case, 
the temperature of the drop is elevated above its normal boiling point. 


Under the usual circumstances when a liquid exposed to a gas phase 
(the atmosphere) is raised to its boiling point it will simply boil and 
maintain a constant temperature. Bubble formation during boiling results 
from nucleation on the walls of the vessel, or from impurities in the boiling 
liquid where microscopic irregularities trap or adsorb gas and present a 
gas-liquid interface where a bubble can be formed. Even this process is 
sometimes not sufficient to prevent superheating, and most chemists are 
familiar with the phenomenon of "bumping" while attempting a distillation. 
The use of boiling chips with rough porous surfaces can prevent this by 
introducing a large gas-liquid interface. Kinetic analysis shows that the 
probable formation of a bubble within a liquid (homogeneously) at the 
normal boiling temperature is very low. Small bubbles tend to collapse 
before they reach a critical size where surface tension effects are not so 
important. 


In the case where an immiscible drop of liquified gas is placed into 
another fluid at a temperature in excess of the boiling point of the liquified 
gas, the drop can remain stable even at remarkably high temperatures. 

For example, several authors have been able to reproduceably heat n pentane 
to 146-147°C at atmospheric pressure (Reid, 1976). This is over 100 echin 
excess of its normal boiling point. 


There is, however, a definite limit to which a liquid can be heated 
called the "superheat" or "cavitation" limit, above which a superheated 
liquid can explode with violence. 


Since the boiling point of chlorine is -34.5°C at atmospheric pressure, 
and the temperature of seawater in Georgia Strait may be expected in the 
5-15°C temperature range, the injection of liquid chlorine into the water 
will result in the formation of superheated drops. It is important to 
ascertain the liklihood of an explosion or explosions following rupture of 
the tank cars. This discussion applies only to that part of the water 
column where chlorine gas is thermodynamically stable. At depths where 
chlorine would normally remain a liquid, gaseous explosion will never be a 
problem. 


Thermodynamic stability analysis shows that for a pure liquid the super- 
heat limit occurs when (Beegle et al-, 1974) 
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9P 
Gy p70 (33) 


Using this thermodynamic criterion and the Redlich-Kwong (1949) 
equation of state, Reid (1976) has provided a generalized curve of reduced 
pressure (P/P.) as a function of the reduced superheat limit (To) 

re 
where P,, and T, are the critical pressure and temperatures respectively 
This curve gives a very adequate representation of experimental observations 
on many gases. 


Based on this curve and the critical pressure (76.8 bars) and critical 
temperature (144°C) of chlorine, the water temperature would have to be in 
excess of 100°C before explosion could take place at atmospheric pressure. 
Clearly, this rules out gaseous explosion as a result of mixing chlorine 
with Georgia Strait water. 


Chlorine Release in Seawater 


On August 20, 1976, some small-scale chlorine release experiments were 
performed in Saanich Inlet at a series of depths from 25 m to 145 m. A small 
size cylinder (10 cm diameter x 36 cm length, Canadian Liquid Air) was 
half filled with liquid chlorine and padded to 35 bars pressure with 
nitrogen. The outlet was connected through a short length (1 meter) of 
-3 cm OD stainless steel tubing, acting as a flow restrictor, to a Whitey 
ball valve. After the ball valve, a short piece of .6 cmstainless steel 
tubing was connected with a one way check valve at the end, and a short 
piece of .6 cmOD Teflon was connected last. The ball valve was opened and 
closed outside the Pisces by a Whitey valve actuator with spring return. 
The actuator was connected through a hydraulic reduction valve (14 bars) 
to the Pisces hydraulic system (69 bars) which was operated from within the 
Pisces. Figure 14 gives a schematic diagram of the equipment, and Figure 15 
shows how the system was arranged on the outside of the Pisces. 


During the release experiments visual examination, colour 35 mm 
slides and super 8 movies were used to document the behaviour of chlorine 
at various depths when released into the water. In the following descrip- 
tion it will be assumed that every 10 meters depth of seawater is approxi- 
mately equivalent to 1 bar of pressure so that, for example, 80 meters 
would have a total pressure of 9 bars, atmospheric pressure being included. 
This should not be in error by more than about 2%. Accurate interconver- 
sion of pressure and depth can be made provided salinity and temperature 
data are available (Saunders and Fofonoff, 1976). 


Initially on the dive, the Pisces descended to the bottom of Saanich 
Inlet about 190 m. Although the order for the cylinder specified padding 
the void space with 35 bars (500 psi) of N,, chlorine did not come out at 
this depth. It is possible that the by inde was not padded to sufficient 
pressure or, more likely, some of the No dissolved in the chlorine causing 
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a reduction in pressure according to Henry's law. The first chlorine 
emerged at 145 m (15.7 bars) as green liquid drops, some running down the 
tube as a thick syrup-like fluid and some squirting out and then quickly 
descending in drops about 1/2 cm in diameter (Fig. 16). No hydrate was 
seen although the chlorine rapidly disappeared from sight below the Pisces 
windows before mixing to any degree with the water. As suggested in the 
section on kinetics, the hydrate formation process is transport controlled 
at the chlorine-water interface and slow mixing here probably slowed down 
hydrate formation even though the hydrate is thermodynamically stable 

(Fis 50, 


At 90 meters the chlorine emerged more vigorously (due to reduced 
hydrostatic pressure) and could be observed for a longer time period than 
before. Figure 17 shows the green liquid coming up out of the tube, and 
then rapidly descending. Some yellowish hydrate could be seen forming on a 
few of the liquid chlorine drops, and fine white flakes of hydrate were 
clearly seen, dispersed through the water a short time after the chlorine 
had been released. These flakes appeared to have almost neutral buoyancy 
neither rising nor falling at an observable rate. The density of the hydrate 
is 1.23 g/cc which indicates that it could eventually sink at a rate depen- 
dent on particle size. Figure 17, also taken at 90 m, shows more clearly 
the liquid chlorine-hydrate clumps which formed shortly after injection 
into the water. These clumps slowly sunk below the Pisces windows. 


At 60 meters depth the chlorine shot out in a fairly dispersed plume 
and formed hydrate more readily than at 90 meters probably due to turbulent 
mixing, and a finer dispersion of chlorine into small drops. Figure eS] 
shows the chlorine plume, and some of the liquid drops with hydrate attached 
are visible to the right of the tube and level with it. Some small bubbles 
were also seen, but these were probably due to N, being released from the 
chlorine similar to the release of co, on opening a carbonated drink. 


At 45 meters a similar behaviour was observed but the liquid chlorine- 
hydrate clumps. which formed initially had occluded or attached bubbles 
which caused them to rise quite rapidly. At this depth, liquid chlorine, 
hydrate and gaseous chlorine were all seen although the gaseous chlorine 
formed slowly. Fine flakes of hydrate were seen shortly after chlorine 
injection as shown in Figure 20. 


At 25 meters no liquid was observed at all, and the chlorine came out 
as a steady gas plume (Fig. 21). Outgassing and expansion of dissolved N 
may have contributed to this, but the general colour showed that a large 
portion of the gas was chlorine. Since liquid chlorine still remained in 
the cylinder (which was inverted) at the end of the experiment, the gas 
seen at this level was generated from liquid chlorine, and not from gas 
coming directly from the cylinder. No hydrate was seen in the water near 
the injection point, but it could possibly have formed further up in the 
water column as the gas plume rose. How far the gas bubbles rose was 
impossible to tell from inside the Pisces Some gas was sighted near the 
surface, but only when the Pisces had cleared its ballast tanks. There 
was no direct evidence of chlorine at the surface. In accordance with the 
theory given in the section on gaseous explosion, no violent explosions 
were seen at any point in the water column during the chlorine releases. 
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Some experiments on chlorine-seawater interaction were performed by 
Macmillan-Bloedel (R.W. Cunnington, 1975) in which chlorine was released 
into seawater in the laboratory, and into the surface water in Burrard Inlet. 
They used a 25 litre chlorine cylinder inverted, so that only liquid chlorine 
was released (The tank was not pressurized with air or N,.) Gaseous chlorine 
and the hydrate formed during all of the experiments, some of the gas 
apparently reaching the surface from a depth of 20-30 cm. Much of the 
chlorine was tied up in "golf ball"sized lumps of hydrate which came to the 
surface and dissociated into the aqueous phase within a couple of minutes. 
(The surface temperature of the water was probably less than 10°C.) 


The physical behaviour of chlorine in seawater conforms quite closely 
with thermodynamic prediction (Fig. 2-5). Gas forms above 45 meters and the 
hydrate clumps rise due to attached bubbles above the depth. In the event 
of a large spill of chlorine at less than 45 meters, it is impossible to 
give an exact estimate of the behaviour, but the majority of chlorine should 
be tied up as a hydrate before reaching the surface. At very shallow depths, 
some chlorine gas could break the surface and escape to the atmosphere, but 
how much would escape depends directly on the magnitude of the leak, the 
depth of the water and the temperature. Lumps of hydrate buoyed by 
occluded gas could also float to the surface, but would tend to dissolve 
into the water rather than release chlorine to the atmosphere. 


Chlorine Release in the Atmosphere 


During an Ocean Chemistry Cruise on the CSS Vector we had the oppor- 
tunity of observing first hand an accidental release of a mixture and 
hydrochloric acid from FMC Co. at Squamish on November 29, 1976. The Vector, 
which at the time of release was on station about 500 meters off shore 
directly in front of FMC was forced off station for about two hours while 
the chlorine dispersed. It was estimated that a mixture consisting of 
76 kg of chlorine and 4.5 kg of hydrochloric acid was lost during the dis- 
charge. The amount of chlorine lost corresponds to less than 0.1% of that 
contained in one of the tank cars. Winds at the time were light (about 
3 knots) and from the west. Initially, a yellowish-white cloud was seen 
at the FMC site (Fig. 22) but soon it spread across the water toward 
Darrell Bay as an elongated whitish mist (Fig. 23). The whitish mist, or 
remants of it were visible along the eastern side of Howe Sound for over 
three hours after the release (Fig. 24), and on returning to station after 
two and a half hours there was still a strong chlorine smell in the air. 
Since the chlorine which was released was mixed with some hydrochloric acid, 
the whitish mist was most likely caused by the interaction of HCl with 
water vapour in the air. Pure chlorine when released might not be so 
visible, although water vapour could lead to the production of HCl by way of 
Equation 3. 


Two researchers on a field trip were conducting a sampling program 
throughout the night (November 29-30) on the tidal flats in front of FMC. 
They both detected chlorine in the air at times, and reported that on 
occasions it smelled strongly while at other times they did not notice it 
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depending on the wind direction. Both developed sore throats as a result 
of the exposure (Thomas and Erickson, 1976). 


Toxicity (Marine Environment) 


As has been shown in the previous sections, the chemistry of chlorine 
in seawater can be fairly complex, and prediction of reaction products in 
the natural environment is virtually impossible (Mattice and Zittel, 1976). 
In spite of this, it appears that toxicity can be fairly well predicted 
provided the correct type of information is available. 


After water has been chlorinated, only three basic forms of chlorine 
need be considered for toxic effects. 


1. Free residual chlorine existing as either HOCI, OCl or even as 
free Cl. The equilibrium between these three species has already 
been dealt with in some detail, and it suffices here to repeat 
that the system is very pH dependent (Fig. 9). 


2. Combined residual chlorine, specifically the chloramines but 
also compounds formed with more complex nitrogenous substances 
such as proteins and amino acids. 


3. Apparent chlorine demand chlorine. This category covers that 
chlorine which has been removed from the system by redox reactions 
resulting largely in the formation of chlorides. Operationally, 
it is the difference between the chlorine added, and the "active" 
chlorine determined by some suitable technique. 


Of these three forms the third contributes almost nothing to toxicity 
in seawater. In the case of polluted waters containing high concentrations 
of anthropogenic materials some toxic secondary products may be formed, 
for example cyanogen chloride (Allen, et al., 1946, 1948) chlorophenols 
(Parker, 1935; Hopkins and Bean, 1966) 5-chlorouracil and 4-chlororesorcinol 
(Jolley, 1973). These compounds should not be important in open ocean water 
and it will be considered here that chlorine removed by chlorine demand has 
been effectively destroyed. 


A further consideration peculiar to seawater is the possibility of 
producing brominated compounds due to the displacement of bromide from the 
seawater by chlorine. Toxicity studies of chlorine in seawater have tradi- 
tionally not been concerned about the generation of bromine, but operation- 
ally this should not cause problems since toxicity measured as a function 
of chlorine addition will cover combined effects of chlorine and bromine. 
The bromide in seawater may indeed explain observed differences in toxicity 
of chlorinated fresh water and chlorinated seawater. In the remainder of 
this section reference will be made only to chlorine while it is realized 
that some or all of the chlorine may be converted to bromine in seawater. 


The combined and free residual chlorine cause almost all of the toxic 
effects in natural waters so that it is important to consider chlorine 
remaining in both of these two forms after chlorination, rather than the 
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original dose of chlorine added (Merkens, 1958; Doudoroff and Katz, 1950). 


Whether the chlorine is in the free or combined form does not seem 
to matter greatly as far as toxicity is concerned, with both free chlorine 
and the chloramines exhibiting the same order of magnitude of effects 
(Brungs, 1973). Evidence has been presented that free chlorine is more 
toxic (Merkens, 1958; Rosenberger, 1971) and that the chloramines are more 
toxic (Holland, et al., 1960). The experiments were not strictly comparable, 
and it is probable that other environmental factors must be considered in 
determing which form is more toxic. The concentration of total residual 
chlorine (free and combined) appears to be both necessary and sufficient 
to determine toxicity (Brungs, 1973; Mattice and Zittel, 1976). (The total 
residual chlorine is best determined by the amperometric, or ferrous DPD 
method. The orthotolidene method has been critized as inaccurate by 
Mattice and Zittel, 1976.) 


So far we have mentioned that the total concentration of available 
chlorine is most important when considering toxic effects. Of equal 
importance is the time of exposure. Therefore, when talking about toxicity 
of chlorine, it is most important to consider the dose-time anticipated. 
Mattice and Zittel (1976) in their recent review have taken all of the 
literature data and have constructed a most useful dose-time curve for the 
toxicity of chlorine both in fresh water and in seawater. Acute and chronic 
toxicity thresholds have been estimated from a statistical analysis of the 
compiled data, and are shown in Figure 24. The chronic toxicity threshold 
is the level of chlorine below which no toxic effects result regardless 
of exposure time. The acute toxicity threshold is constructed from considera- 
tion of dose-time effects, where high concentrations can be endured for 
shorter time periods. There is probably a limit on the acute threshold 
curve where higher concentrations of chlorine do not shorten death time. 


From Figure 24 it can be seen that the chronic toxicity threshold 
of chlorine in the marine environment is about .02 mg/l. Since the thres- 
hold is a statistical creation, we cannot say that there will be no toxic 
effects from concentrations below this, but we can say that the chance of 
mortality is very small, and becomes smaller the further below the thres- 
hold line the total residual chlorine concentrate is. 


Exposure to higher chlorine concentrations can be endured by a wide 
spectrum of life forms for a short time period with a .1 mg/l level for 
one minute not contributing too great a hazard. 


These toxicity thresholds have been generated by reference to tests 
on many different animals in a variety of developmental stages. As a result 
the thresholds tend to represent the safety limit for the most sensitive 
species, and will in general err on the side of safety. 


Toxicity (Atmospheric) 
Chlorine is an extremely irritating substance to mucous membranes and 


the respiratory system. The minimum concentration causing slight irritation 
is about 1 1/1 (Matheson Gas Data Book, 1966) but chlorine cannot be detected 
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in the atmosphere by smell until it reaches a level of about 3.5 ul/1. 

The maximum level which can be breathed for one hour without causing damage 
is 4 1/1. A detectable odour of chlorine therefore indicates that there 

is a safety hazard, and steps should be taken to reduce chlorine inhalation. 
At a level of 15 pl1/1 there is an immediate throat irritation, and levels of 
50 ul/1 are dangerous even for short exposures (less than 30 minutes) 
(Dangerous Properties of Industrial Materials, 1975). Levels of 1000 yl/1 
are likely to be fatal after a few deep breaths. Chlorine when inhaled 
reacts with water to form HCl and nascent oxygen both of which attack tissue 
and can result in pulmonary edema. 


Conclusions 


1. Thermodynamics show that four phases may be expected when chlorine and 
water are mixed, gaseous chlorine, liquid chlorine, chlorine hydrate 
(solid) and water. With sufficient dilution chlorine will dissolve 
completely in water forming a "bleach"'. Above 50 meters gaseous 
chlorine can form at the anticipated range of temperatures. 


2. Chlorine interacts with water on dissolving to form several species. 
Cl., HOC1, and OCl . The equilibrium involving these species will be 
important in considering the effect of chlorine on pH, alkalinity and 
even toxicity. 


3. Chlorine solutions in seawater will be much more strongly oxidizing than 
normal oxic seawater, and it is likely that the 67.3 mg/1 bromide in 
35°/oo salinity seawater will be converted to bromine. 


4. Most of the chlorine speciation reactions are rapid, and hydrate forma- 
tion is very fast,being transport controlled. It is, therefore, 
anticipated that given mixing, most of the chlorine will quickly form 
hydrate which can then dissolve on further dilution with water. Even if 
tank car rupture occurs in shallow water, most of the chlorine should 
be tied up as a hydrate. 


5. The chlorine demand of seawater is about 1.5 mg/l (half hour) and 
about 3 mg/l after two days. Dilution of chlorine to 3 mg/1 should 
result in it being removed via redox processes within a couple of days. 


6. Chlorine release experiments show that the chlorine-seawater system 
conforms very well to thermodynamic prediction, at least for small 
releases. 


7. Gaseous explosion caused by localized superheating of liquid chlorine 
is not expected to be a problem at the anticipated water temperatures. 


a0 


Chlorine is expected to exhibit toxic effects when present above 
.02 mg/1, but should not cause problems below that level. Once 
the chlorine has been removed by the natural demand of seawater 
it is effectively out of the system. Secondary toxic products 

should not become a problem in a relatively unpolluted area such 
as Malaspina Strait. If the tank cars are in deep water, leaked 
chlorine should settle to the bottom and toxic water formed from 
mixing the chlorine with seawater should thereafter remain in 


deep water until sufficient dilution removes the chlorine below 
toxic levels. 
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Conversion of Common Units to SI Units 


To Convert From Common Unit 
Atmosphere (Atm) 

Calories (Cal) 

pounds/sq inch (psi) 


Tons (long) 
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TABLE I 


Toesl, UnLe 
Bar 
Joules (j) 
Bar 


Metric Tons 


Multiply By 
13:02:35 
4,184 

06895 


i LOG 


TABLE II 


Selected Properties of the Chlorine Water System 


Chlorine 

Molecular Weight 

Boiling Point 

Specific Gravity (liquid) 0-20°C 
Latent Heat of Vaporization 28.3°C 
Specific Heat (liquid) 


Specific Heat Cp (gas) 


Hydrate 


Composition (Cl, - n H,0) 


2 
Critical Decomposition Temperature 


Specific Gravity 


Heats of Reaction 


a) cl, Ce eran H,0 (lL) > cl, 


70.9 g/mol 
=34755-C 
1.47-1.41 
60.6 cal/g 
0.226 cal/g 


0.114 cal/g 


“on H,0 (Cal/mol at*28273°C) 16,550 


b) Cl, (4) +n H,0 Do” Cl, - on H,0 (Cali/molaat 2873336 %022 9253 


c) Cl, (LZ) + H,0 (Ll) > HC10. 5 + HCl. (Gal/moltat 25°¢):a74 5420 


IMatheson Gas Data Book 
2Fernandez, et al. 


3CRC Handbook of Chemistry and Physics, 50th edition, 1970. 
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Malaspina Strait (see Fig. 1). 
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seawater, and the redox potential as a function of chlorine added. 
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with an initial pH of 7.9, and at total carbonate concentration from 
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Schematic diagram of the equipment used to release small amounts of 
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A picture of the equipment on the outside of the Pisces. 
Liquid chlorine drops sinking in seawater at about 145 meters depth. 


Liquid chlorine drops sinking in seawater at about 90 meters. Yellowish 
hydrate can be seen forming on some of the drops. 


Chlorine in seawater at 90 meters a short time after injection. 
Clumps of liquid chlorine and hydrate can be seen in the water. 
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Chlorine plume at 60 meters. Liquid chlorine drops with hydrate 
attached are visible to the right of the tube above and below the 
nozzle. 

Fine flakes of hydrate generated shortly after chlorine release at 
45 meters. Flakes similar to these were also seen at 60 meters and 
90 meters. 


Plume of gaseous chlorine at 25 meters. 


Chlorine "cloud" at Squamish shortly after an accidental discharge 
from FMC. 


Elongated chlorine cloud about half an hour after the release. 


Chlorine mist hanging above the waters of Howe Sound about three 
hours after the release. 


Acute and chronic toxicity thresholds for chlorine in seawater. 
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Figure 5. Phase diagram for the chlorine-water system at the hydrate com- 
position. 
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Solubility of chlorine in seawater as a function of pressure and 
temperature. 
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Figure 8. In situ pH measured as a function of depth for two stations in 
Malaspina Strait (see Fig. 1). 
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Figure 23. Elongated chlorine cloud about half an hour 
after the release. 


Figure 24. Chlorine mist hanging above the waters of Howe 
Sound about three hours after the release. 
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ABSTRACT 


Physical, chemical and biological oceanographic observations are made from 
the weathership at Ocean Weather Station Papa, and between Esquimalt and 
Station Papa, on a routine continuing basis. Physical oceanography data only 
are shown, including surface observations and profiles obtained with bottle 
casts and conductivity-temperature-pressure instruments. 
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INTRODUCTION 


Canadian operation of Ocean Weather Station P (Latitude 50°00'N, 
Longitude 145°00'W) was inaugurated in December, 1950. The station is occu- 
pied primarily to make meteorological observations of the surface and upper 
air and to provide an air-sea rescue service. The station is manned by two 
vessels operated by the Marine Services Branch of the Ministry of Transport. 
They are the CCGS Vancouver and the CCCS Quadra. Each ship remains on station 
for a period of six weeks, and is then relieved by the alternate ship, thus 
maintaining a continuous watch. 


Bathythermograph observations have been made at Station P since 
July 1952. A program of more extensive oceanographic observations commenced 
in August 1956. This was extended in April 1959, by the addition of a series 
of oceanographic stations along the route to and from Station P and Swiftsure 
Bank. These stations are known as Line P stations. The number of stations on 
Line P has been increased twice and now consists of twelve stations (Fig. 1). 
Bathythermograph observations and surface salinity sample collections, in 
addition to being made on Line P oceanographic stations, are also made at odd 
meridians at 40', i.e. 139940'W, 141940'W, etc. These stations are known as 
Line P BT stations. Data observed prior to 1968 has been indexed by Collins 
et al (1969). 


The present record includes hydrographic, continuously sampled STP . 
and surface salinity and temperature data collected from the CCGS Quadra 
during the period 3 December 1976 to 13 January 1977. 


All physical oceanographic data have been stored by the Canadian 
Oceanographic Data Centre (CODC), 615 Booth Street, Ottawa, Ontario, Canada. 
Requests for these data should be directed to CODC. 


Biological and productivity data are published in the Manuscript 
Report series of the Fisheries Research Board of Canada (FRB), Pacific 
Biological Station, Nanaimo, British Columbia, Canada. Requests for these 
data should be directed to FRB. 


Marine geochemical data are for the Ocean Chemistry Group, Ocean 
and Aquatic Sciences, Environment Canada, 512-1230 Government Street, Victoria, 
British Columbia, Canada. 


PROGRAM OF OBSERVATION FROM CCGS QUADRA, 3 DECEMBER 1976 - 13 JANUARY 1977 
P-76-9) (CODC Ref. No. 15-76-009) | 


Oceanographic observations were made by Mr. B. Whitehouse of Seakem 
Oceanography Ltd., Victoria, B.C. 


En route to Station P, Line P Stations 2 to 12 were occupied and an STD 


profile made to near bottom or 1500 metres. One hydrocast to 1500 m was done 
at Station 10. 


Samples for salinity, nitrate, nutrient, alkalinity and total CO, were 
taken from the seawater loop at all whole and half stations except 128. 
Surface bucket salinities were taken at all whole and half stations. Surface 
bucket temperatures were taken at all whole stations. 

Surface tarball tows were made at Stations 2, 4, 6, 8, 10 and 12. 


The thermosalinograph, surface temperature recorder and PCO, system were 
run continuously (though the thermosalinograph produced erratic traces). 


Mechanical BT's or XBT's were taken at all whole and half stations. 


No samples were collected at Station 12% since the bridge did not inform 
the oceanographer of the ship's position. 


At Station P the oceanographic program was carried out as follows: 


I. Physical Oceanography 


1) Profiles of salinity, temperature and oxygen were obtained from 3 hydro- 
casts to near bottom (4200 metres); one cast was to 600 m only. 


2) 6 STD profiles to 1500 metres and 14 to 300 metres were obtained. 


3) BT's or XBT's were taken every three hours to coincide with meteorological 
observations, encoded and transmitted according to the IGOSS format. 


4) Salinity samples were collected daily at 0000 hrs GMT from the seawater 
loop. 


Il. Marine Geochemistry 


1) Nutrient and salinity samples were collected daily at 0000 hrs GMT from 
the seawater loop. One 24 hour series of nutrient samples was taken 
each hour from the seawater loop. One profile for nutrients to 500 m 
and one profile for tritium to 500 m were taken. 


2) 


3) 
4) 


5) 


6) 


7) 


ELL 


a) 


2) 


Alkalinity and total CO) samples were taken every 2 to 4 days from the 
seawater loop and in addition 2 profiles to 500 m were taken. 


Air CO) samples were taken in quadruplicate at weekly intervals. 
4 surface tarball tows were completed. 


3 seawater C-14 samples were extracted from 45 gallons of seawater taken 
from the seawater loop along with 3 seawater C-13 and 3 air C-13 samples. 


6 hydrocarbon samples were obtained. 


PCO2 carboys were filled 10 times while on station. 


Biological and Productivity 


Samples were obtained as follows: 


11 - 150 metre vertical plankton hauls. 
1 - 300 metre vertical plankton haul. 
2 - 1500 metre vertical plankton hauls. 


2 profiles to 200 metres for each of plant pigment, nitrate and C 


productivity were obtained. 


En route from Station P, an STD was made at Stations 12 to 1 and one 


hydrocast was done at Station 6. Nutrient, nitrate, alkalinity and total CO, 
samples were taken from the seawater loop at Stations 12 to 1. Salinity 
samples were taken at all whole and half stations 12% to 1. Surface bucket 
temperatures were taken at all whole stations. Tarball tows were taken at 
Stations 12, 10, 8, 6, and 2. Mechanical BT's or XBT's: were taken at all 
whole and half stations. 


Observations for Other Agencies 


1) 


2) 


3) 


Marine mammal observations were made by the ship's officers for Mr. I. 
McAskie, Fisheries Research Board of Canada, Pacific Biological Station, 
Nanaimo, B.C., Canada. 


Bird observations were made by the ship's officers for Mr. M. Myres, 
University of Alberta, Calgary, Alberta, Canada and Mr. J. Guiguet, 
Curator of Birds and Mammals, Provincial Museum, Department of Provincial 
Secretary and Travel Industry, Victoria, British Columbia, Canada. 


Air CO, samples weekly in duplicate for Scripps Institution of Oceano- 
graphy, La Jolla, California, U.S.A. 


Data was processed for publication by Ms. M. Sainsbury of Seakem Oceano- 


graphy Ltd:,,, Victoria, B.C. 


OBSERVATIONAL PROCEDURES 


Observations for salinity, oxygen and temperature from all hydrographic 
casts, including the surface, were obtained with Niskin water sample bottles 
equipped with either Richter and Wiese and/or Yoshino Keiki Co. reversing 
thermometers. Two protected thermometers were used on all bottles and one 
unprotected thermometer was used on each bottle at depths of 300 m or greater. 
The accuracy of protected reversing thermometers is believed to be + 0.02°C. 


The daily surface water temperatures were measured from a bucket sample 
using a deck thermometer of + 0.1°C accuracy. The daily surface salinity 
samples were obtained from the seawater loop. When the seawater loop was 
not operational these samples were obtained with a bucket, and are indicated 
with a 'b' in this data record. 


Salinity determinations were made aboard ship with either an Autolab 
Model 601 Mark III inductive salinometer or a Hytech Model 6220 lab salino- 
meter. Accuracy using duplicate determinations is estimated to be + 0.003 /oo. 


Depth determinations were made using the ''depth difference" method 
described in the U. S. N. Hydrographic Office Publication No. 607 (1955). 
Depth estimates have an approximate accuracy of + 5 m for depths less than 
1000 m, and + 0.5% of depth for depths greater than 1000 m. 


The dissolved oxygen analyses were done in the shipboard laboratory by a 
modified Winkler method (Carpenter, 1955). 


Line P engine intake continuous temperature on both ships were recorded 
by a Honeywell Electronik 15 Recorder. The temperature probe is at a depth 
of approximately 3 metres below the sea surface and the instrument accuracy 
is believed to be + 0.1 C. 


Each ship is equipped with a Plessey Model 6600-T thermosalinograph which 
is used, on Line P, for continuous recording of surface temperatures and 
salinities from ship's seawater loop. The temperature probe is mounted at the 
seawater loop intake (approximately 3 metres below the surface) and the 
salinity probe and recorder are situated in the dry lab. The accuracy of this 
instrument is believed to be + Q.1°C for temperature and + 0.1°/00 for salini- 
ty. 


STD profiles were taken with a Plessey Model 9006 STD system. 


COMPUTATIONS 


All hydrographic data were processed with the aid of an IBM 370 computer. 
Reversing thermometer temperature corrections, thermometric depth calculations 
and accepted depth from the "depth difference" method were computed. Extran- 
eous thermometric depths caused by thermometer malfunctions were automatically 
edited and replaced. A Calcomp 565 Offline Plotter was used to plot temper- 
ature-salinity, and temperature-oxygen diagrams, as well as plots of tempera- 
ture, salinity, and dissolved oxygen vs log depth. These plots were used to 

10 
check the data for errors. 


Missing hydrographic data were obtained using a weighted parabolas inter- 
polation method (Reiniger and Ross, 1968). These data are indicated with an 
asterisk in this data record. 


Data values which we suspect but which we have included in this data 
record are indicated with a plus. These data have been removed from punch 
card and magnetic tape records. 


Analog records from the salinity-temperature-pressure instrument have been 
machine digitized, then replotted using the Calcomp plotter. 


Digitization was continued until original and computer plotted traces 
were coincident. Temperature and salinity values were listed at standard 
pressures; integrals (depths, geopotential anomaly, and potential energy 
anomaly) were computed from the entire array of digitized data. 


The headings for the data listings are explained as follows: 


PRESS is pressure (decibars) 

TEMP is temperature (degrees Celsius) 

SAL is salinity (parts per thousand) 

DEPTH is reported in metres 

SIGMA-T is specific gravity anomaly 

SVA is specific volume anomaly 

THETA is potential temperature (degrees Celsius) 

SVA (DEPTH) is potential specific volume anomaly 

DELTA D is geopotential anomaly (J/kg) 9 

POT EN is potential energy in units of 10 ergs/cm 

OXY is the concentration of dissolved oxygen expressed in millilitres 
per litre 

SOUND is the velocity of sound in m/sec 
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Figure 2. Composite plot of temperature vs 1lo8ig depth for Line P 
stations. P-76-9. 
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Figure 3. Composite plot of salinity vs lo8io depth for Line P 
stations. ~ P-/6-9. 
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Figure 4. Composite plot of temperature vs 10819 depth for Station P. 
P-/6-9. 
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Figure 5. Composite plot of salinity vs 10819 depth for Station P. 
P-76-9. 
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Figure 6. Composite plot of oxygen vs 10819 depth for Station P. 
P-/6-9. 
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REFERENCE NUe 76 9= 17 UATE Bees tes 7c STATIUN P 
PCSITICN SO DOeONe 145= 0 -0W GMT 1763 
RESULTS CF STP CAST 6S POINTS TAKEN FRUM ANALLS TRACE 
PRESS TEMP SAL C&ePTH SIGMA SVA Della POT. STUINO 
yj v oa 
0) Oe SS 32264 @) 2560S pede ho yr yl ©: Oyen Oe0 1474. 
10 6052 32264 10 253605 235 +3 e224 Oo. OF LWT 4. 
20 Geoo 32-64 20 25005 22535 Ce 47 OS 1474.6 
30 5253 32264 30 256 OS ZS eG OeZ7l Oelil l Aet Sie 
50 Ge Os 322004 510 ct365 23 res 1els Ote3 © ) tery Bis 
75D Se27 32275 75 255/59 Pisa e Lb o¥& Oe om | | 2a? Li's 
100 4e21 33203 99 2boe2ea 181.3 2ec4 Leal POW ec 
| at | 4252 33254 124 €€e50 14623 Jie 6 Le Seé 14O0G 6 
150 4255 33072 149 2oe74 Laas 3 Se00 Z2iewis L (i7Ore 
175 4e34 33-79 174 26281 letee Smee nee 1 2a? Cu 
200 4e18 33281 199 26234 123660 ey 9). Sroul es 146096 
Va ds 4207 33284 2cu 26038 120.0 3eS4% 30 84 14096 
Bes 6) 3299 aS YG he Ke Wy 6 248 26e%1 PbiwredD 4023 405% L4#oGe 
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Hi 
Ch 
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SALINITY, 0/00 


OFFSHORE OCEANOGRAPHY GkKOQUP 


REFERENCE NOe 


PCSITICN 


PRESS 


762 Je 


50—= OeON>e 
RESULTS OF STP CAST 


Temp 


6058 
6258 
6258 
6258 
6e58 
0e58 
4022 
4256 
4260 
4226 
4214 
4207 
4.200 


SAL 


32262 
32e¢62 
32262 
32263 
32663 
32265 
3291 
33244 
Set 2 
33076 
33-79 
33-82 
33-86 
33-91 


19 
145— 0.0W 


CEPTH 


¢) 
10 
20 
30 
50 
75S 
99 
124 
149 
174 
199 
223 
248 


61 


DATE 159/712/77€ 

GMT 1723 

71 POINTS TAKEN FROM 

SIGMA SVA 

T 
250663 237 22 
25263 237¢5S 
25e62 Zhe 
252663 237-40 
22¢OS 232 
2Se65F 23E-0 
26013 190 64 
2€e5l 154.3 
26073 134.0 
26077 i 30-0 
26033 12563 
26286 L22063 
26290 1186S 
2060094 114.7 


298 


ANALUCG 


DELTA 
v 
Cel 
Oec4 
Co47 
Oe?l 
1.219 
1.78 
2030 
eel 
3eC8 
3041 
see 3 
4004 
4ei4 
49392 


STATIGN P 

TRACE 

POT e SIIIND 
Ged 1474. 
OcOLt  te7s. 
Oe05 1475. 
Ovlf L£eF7Se 
O«30 2475S 
Ovee -TH766 
Reta 467. 
1663 140596 
2et3 14704 
Ze67 TS7igk 
Jel. £E7Ce 
69S £670. 
40268 1470. 
SeIh tare. 
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DB 
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C 
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TEMPERATURE, € 
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MO.-12 DAY-16 GMT-17.3 
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OFFSHORE OCEANCGRAPFY GROUP 


REFECRENC 
PCSITION 


c NOe6 
5 0= 


76= 9J= 
OeON, 


RESULTS GF STP CAST 


PRESS 


0 

10 
20 
30 
SO 
rf.) 
100 
izs 
150 
175 
200 
22S 
250 
300 
400 
S00 
600 
800 
1000 
1200 


TEMP 


6057 
6057 
6057 
6057 
6057 
6256 
4238 
4.53 
4268 
9eS52 
4231 
4218 
4208 
4203 
3e 36 
3269 
3250 
3013 
2082 
2256 


SAL 


32e¢62 
32263 
32263 
32264 
3263 
32264 
32-90 
33244 
33-71 
33-76 
33-80 
33482 
332484 
33491 
34204 
34e13 
34220 
34e31 
34238 
346 44 


2l 
145— 0O.0w 


101 PUINTS TAKEN FROM ANALUG 


DEPTH 


O 
10 
20 
30 
50 
7S 
99 

124 
149 
174 
199 
223 
248 
298 
337 
496 
995 
793 
990 
1188 


63 


CATE 16/7172/76 


GMT 1723 


SIGMA 
T 
25263 
25-63 
25 264 
25204 
25264 
25004 
26210 
26052 
26-72 
26079 
26e82 
26 e385 
260-88 
26294 
27-06 
27.2158 
27 e223 
27235 
27443 
27edl 


SVA 


237-0 
23608 
2sGeS 
236061 
23701 
23607 
1G2eY9 
1539 
135.3 
128638 
12564 
L228 
12026 
11iSe0 
1044 

S606 

GOeo«l 

7906 

V20¢5 

6507 


DelTa 
Db 
0209 
GCeed 
0247 
Oe ll 
leld 
1.78 
2031 
207d 
Sell 
3045 
Sel DS 
4206 
4237 
4eSe 
OeO05D 
72Go 
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326d 
1120 
heool 


STATiGn P 
TRACE 
PUT. SJIUNO 
iN 
0 60 1474. 
Oecl 1474. 
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Oell 1475. 
O2.22C A47Se 
0263 1475. 
1eid 1457. 
Levee 1140S. 
2eid 1470.6 
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Bie S10 147Ce 
3048 147Ce 
4eo7l1 14706 
G6e3s 147Ce 
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o 
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UFFSRFORE QUCEANCGKAPFY GRCUP 
REFERENCE NGe 
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76=> 9= 
OeONe 


RESULTS CF STP CAST 


PRESS 


QO 
10 
20 
30 
SO 
75 

100 
izs 
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i7s 
200 
225 
2sc 


TEMP 


6eS1 
6e5l 
6e51 
605i 
6e52 
6052 
4219 
4232 
4206 
4254 
4237 
4e22 
G4el2 


SAL 


22663 
32.63 
32263 
32064 
32263 
32264 
32.89 
3B3e1€ 
33265 
33-77 
33.80 
33-81 
33.82 


23 
145— 02.0Ww 


71 PCINTS TAKEN FROM AWALES 


DEPTH 


fe) 
10 
20 
30 
50 
75 
99 
124 
149 
174 
199 
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248 
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DATE 17/12/76 


GMT 17.3 


SIGMA 
T 
25264 
25264 
25265 
25265 
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25865 
26e12 
2€234 
26-67 
26278 
2€e82 
26234 
26237 
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22506 
23509 
23507 
23504 
23662 
22661 
151.3 
i7o.9 
139.9 
12947 
125.9 
123.8 
12104 


DLL TA 
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1474. 
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OFFSHORE OCEANOGRAPKY GRCUP 


REFERENC 
POSITION 
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S0O= 


76> G= 
OeONe 


RESULTS CF STP CAST 


PRESS 


TEMP 


62043 
6243 
6243 
6044 
6044 
Sedl 
4224 
4255 
4273 
4047 
4034 
4-17 
4ell1 


SAL 


32062 
32-62 
3262 
32062 
32262 
32070 
33202 
332523 
33-73 
33077 
33-79 
33280 
33285 
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DATe 18/712/7E 


GMT 172 
69 POINTS TAKEN FRUM 


SIGMA 
T 
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26e5& 
26272 
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OFFSHORE QCEANCGRAPHY GROUP 


KCFERENCE NCe 
50] 


POSITION 


162" Ger 25 
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RESULTS OF STP CAST 


PRESS 


0 
10 
20 
30 
50 
75 

100 
1<cS 
150 
175 
200 
225 
250 
300 


TEMP 


6241 
6041 
6040 
6040 
6041 
6239 
4e27 
439 
4.57 
4241 
4e21 
4213 
4209 
4203 


SAL 


0-20W 
71 POINTS TAKEN FROM ANALUS 


DEPTH 


124 
149 
174 
1Sg9 
225 
248 
298 


69 


CATE 19/12/76 


GMT 17.23 


SIGMA 
T 
25265 
2520S 
€5e65 
<cSe65 
25057 
260613 
26248 
260074 
26280 
26284 
26288 
26091 
26097 


SVA 


23 De1 
23504 
2355 
28505 
2350S 
234-20 
1S0e5 
197 es 
1333 
leé7ee 
1235 
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DELTA 
v 

Oe 
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3041 
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TRACE. 
PUT. SUUND 
EN 
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71 
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STATLON P 

TRACS 

POT e SIUND 
éN 
Oe VD 1472.6 
OeVl 1373. 
0205 1474. 
Vell 1474. 
Oe 30 1474. 
Oc? 1474. 
bie k2 14607. 
1etO 140Ge 
2007 1470. 
2esG 1495S e 
3ei1¥9 140%3-6 
Seis 1453.6 
4ede 147Ce 
Oecd0 1470.4 
10e0S Ll472-656 
l4e7e 1473. 
1S eSG 1473-6 
ERG (Nd 1475. 
46e2e 1478. 
&2ee0G 148306 


DB 


Maer tigiee 


60 


120 


—— 
OO 
&) 


20: 


300, 


Bie 


PEMP ERA RUSE 
S 8 


REF. NO. 76 - > 


DO-O.O0 N 145-0.0 W 


MO.-12 DAY-22 GMT-17 


oo 
1 


SAL INI 


G 
Le 


©: 
, O70 


16 


3. 


73 


OFFSHORE QCEANOGGRAPHY GRGUP 


REFERENCE NOe 76= 9= 28 DATE 22/12/76 STATICN P 
PUSITION SO= OeON, 145=— 0.0W GMT 17.23 
RESULTS GF STP CAST 63 POINTS TAKEN FRUM ANALUS TRACE 
PRESS TEMP SAL DEPTH SIGMA SVA DOi.LTA PUT. SUUND 
T D miN 
0) 6220 32<63 0 25068} | 2316S Oe oerne) 14736 
10 6220 3263 10 25268 eee Coe O«eCil 1473. 
20 6e20 32263 20 25.68 2322.3 0246 Oe0S 1473. 
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| 4248 33057) 124 26063 143.3 2eG3 Liewsé 1454S. 
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OFFSHORE OQCEANCGRAPHY GROUP 


REFERENCE NOe 


PCSITICN 


100 
leo 
150 
i75 
200 
Fy ies 
250 
300 


76> G= 


50= Oe0ONe 
RESULTS OF STP CAST 


TEMP 


6¢35 
6e¢35 
6¢e35 
6e35 
6633 
6027 
4216 
4239 
4252 
426 
4220 
4210 
4206 
4200 


SAL 


32264 
32264 
32264 
32.264 
322584 
32067 
33202 
3351 
33-75 
33280 
33-85 
33287 
33-89 
33-96 


29 
145—= 0-0W 
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CEPTH 
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io 
20 
30 
So 
75 
99 

124 
149 
174 
1Ss9 
223 
248 
298 


75 


DATt 25/12/7€ 


GMT 174 


SIGMA 
1 
25067 
25007 
e567 
25267 
25-67 
25070 
260622 
2E€e¢59 
26276 
2£€283 
26287 
26290 
26492 
26098 


SVA 


2328 
essec 
egze3 
23344 
233065 
23028 
1381-5 
147.23 
V3lei 
124e7 
120¢5 
118-1 
11662 
1106G 


DELTA 
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Oe 
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Oe7VU 
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1.75 
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STATIOGN P 

TRAGS 
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3819 1409. 
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KEFERENCE NOe 


PCSITIGN 


PRESS 


76> Go 


49—284e0N> 
RESULTS OF STP CAST 


TEMP 


6e16 
6216 
6216 
6el7 
6e15 
4260 
4265 
4253 
4033 
4219 
4-07 
4-200 
3295 
3e81 
3066 
3e49 
3215 
2e82 
2059 


SAL 


32264 
32264 
32264 
32264 
32-65 
32265 
33204 
33052 
33-75 
33079 
33-82 
33285 
33.87 
33293 
34204 
34214 
34e21 
34231 
3440 
34244 


30 
145—48 .0Ow 


CEPTH 


) 
10 
20 
30 
50 
75 
99 
124 
149 
174 
199 
223 
248 
298 
397 
496 
595 
793 
990 
1188 


77 


CATE 28/12/7€ 


GMT 17.6 
122 POINTS TAKEN FROM ANALCG TRACS 


SIGMA 
T 
<5e70 
25270 
25470 
2570 
25270 
25e7C 
2601S 
2€ 056 
26077 
26032 
2€.85 
2668S 
26091 
26097 
Bw Oit 
27016 
27224 
27434 
27044 
27630 


SVA 


23006 
230¢eS 
231.20 
231iel 
2420649 
23C eS 
184.3 
1494.5 
130-5 
126.61 
12267 
119.3 
117.6 
1S2e7 
G5e6 
8Sel 
7928 
71lel 
6020 


CLL T & 
1D) 
Oe’ 
Vels 
Oe 4H 
0609 
1215 
Wee eS 
Z2eedD 
2eF7 
3201 
39 35 
3e€4 
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4e8e 
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Oe 8 
VeEe 
9e4y 
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22. 37 


STATIGN & 


PT e 
eN 
00 
0.-O0l 
0205 
Qeli 
Vecs 
Oe EE 
lele 
1.260 
2eV0E 
PAL eNs | 
See l 
32d6 
4e5F 
Hels 

10.C5 

14.59 

1S.77 

3iebG 

45448 

BiOrenGe 


S IUND 


1473-6 
14736 
PAW ioe 
1473. 
1473 6 
1474. 
146383. 
147C. 
147C. 
1460S. 
1460S. 
i4UGe 
149G.6 
147C. 
1471. 
1472. 
1473. 
14756 
1477.6 
1480. 


DB 
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G3 


MN 
ne 
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90-0.0 N 145-0.0 W 
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TT 
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S 
, €700 


16 


52 


OFFSHORE OCEANCGRAPHY GROUP 


REFERENCE NOe 
POSITION 


5'0= 


RESULTS CF STP CAST 


PRESS 


QO 
10 
20 
30 
50 
75 
100 
i2s 
150 
175 
200 
225 
250 
300 


TEMP 


6007 
6207 
6207 
6207 
608 
5Se77 
4e20 
4e31 
4252 
4238 
4e20 
4e10 
3-939 
3093 


SAL 


32266 
32266 
32266 
32266 
32e66 
32e73 
3296 
33e29 
3369 
33676 
33-81 
33283 
33486 
33-91 


76= 9= 3i 


OeON, 145= O0-0w 


DEPTH 


) 
10 
20 
30 
50 
75 
99 
124 
149 
174 
199 
223 
248 
298 


ie) 


CATE 29/12/76 
GMT 172 
60 PSCINTS TAKEN FRUM 


SIGMA 
T 
eSe7le2 
25e7e2 
25-72 
2S5e72 
25072 
25081 
26017 
26042 
eOel2 
2607S 
26284 
26-37 
2E290 
2609S 


SVA 


22 8:10 
22264 
228 e5 
22806 
22869 
22044 
186.2 
1622.7 
13409 
1285 
12345 
l2lee 
118.2 
114.0 


ANALUS 


Le LTA 


STATION 


TRACE 


PCT. 
EN 
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UFFSRFCRE GCEANCGRAPKY GROUP 


REFERENCE NC. 76= 9= 48 GATED 172 17 77 STATION, 
PCSITIUN 49=10.60N, 132=40.0W GMT 1468 
BeSUutTsS CF STP CAST 112 POINTS TAKEN FROM ANALOG TRACE 
PRESS TEMP SAL DEPTH SIGMA SVA LLLTA POT. SIUND 
T D 2 
0 9265 32038 0 2405S 25766 Ved 0-0 1486.6 
10 9265 32.38 10 2409S 298.20 Ue 3v CeC2 1436. 
20 9265 32.38 20 24699 2982 0260 0.06 1486. 
30 9268 32239 30 2469S 259863 Oe&9 0.14 1486. 
50 9073 3241 50 25200 2598.0 1.49 Oe3& 1487. 
75 7201 32-64 7S 25059 242.3 2018 0e32 1477. 
100 6256 32684 99 25080 22169 2e7& 1633 1476. 
N29 620 JS e30 124 260621 18362 Saat 16SZ2 1476. 
150 6ell 33.66 149 26052 15461 3669 2050 1476-6 
175 599 33.82 174 26065 14264 460% 3eil 1476. 
200 5265 33.88 199 2667S USGS 4041 Sut IGS 
225 5.33 33.87 223 26070 Uwe 4074 4045 1475.6 
250 Se 06 33.89 248 26081 12764 54 06 Se2¢ 1474. 
300 476 33.91 298 26686 122.69 De tG Be62 tats. 
400 4030 33.99 397 26697 11361 0&0 big2a (Tea. 
500 3299 34-08 496 27-08 102-6 70S 160.18 1474. 
600 3.78 34615 595 27016 SOI 3095 21679 1475. 
800 3047 34.29 793 27230 C407 10e7t 34.68 DP gh 
1000 3209 34.38 390 27-41 P50 la he s Su 460432 1475. 
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OFFSHORE OCEANCGRAPHY GROUP 


KEFERENCE NUe F62 9= 51 OATE, 124, 1/77 STATIN 
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SURFACE SALINITY AND TEMPERATURE GCESERVAT IONS 
CRUISE REFERENCE NUMBER 76= 9 


DATE/TIME SALINITY TEMP LONGITUDE 
YR MO DY GMT 0/00 weST 
76 a2 oeesOo 322044 125-33 
76 12> (4 140 31526 126= 0 
765429 4345 31.833 126-40 
2G) dee). AtPo7 OG 3236C 127=40 
76522), Bits 326424 128=40 
76 12 4 1455 320327 129-40 
76 12 4 1755 322400 130-40 
7612 “4 2140 32.394 13140 
76 12 «5 30 326325 132-40 
76 12 5S §20 326352 133-40 
76 2) 5 Base 32 0362 134=40 
76 12 5 1105 32.358 135-40 
76 12 5S 140C 320424 136-40 
76 12 § 1745 326511 137-40 
76 12 5 2040 32-500 138-40 
76 12 "6 30 32.479 139-40 
76 12 6 420 32.523 140-40 
76 12 6 B10 322554 141-40 
76 12 6 1315 32657 Tea 142-40 
76 12 8 0 32 635 6 06 ON STATION 
76 ¥Ye.9 © 3266323 666 ON STATION 
76.12.10 ) 32-626 609 ON STATIGN 
76 121 fe) 326E€4C 604 ON STATION 
76 12 12 ) 32-619 66 ON STATION 
76 12 13 fe) 32 6637 6 04 ON STATION 
76 12 14 ) 32663¢ 606 ON STATION 
76 12°15 O 32.631 665 ON STATION 
76 12 16 ry) 32-2634 605 CN STATION 
T6021 7, fe) 32-643 604 ON STATION 
76. P2418 f) 32 2634 604 ON STATIOGN 
76.12 19 C 3262S 603 ON STATION 
76 ee ) 322632 Ee ON STATION 
7612621 0 32.62 604 ON STATION 
15 iene ) 322634 664 ON STATION 
7h Fel 23 ) 32 2631 Gee ON STATION 
TH 134 C 326625 664 ON STATIGN 
76 12 25 ) 32.605 Set ON STATION 
76, 12 26 0) 32463C Geil ON STATION 
76 42 27 fe) B276a3 006 ON STATIUN 
76)12e28 fe) 32 6637 605 GON STATIUN 
76° 12 39 ) 326635 620 GN STATION 
76°12 Bo ) 32.639 620 ON STATION 
TEAS St re) 326638 60 CN STATION 
af ene? Calon | ) 32.581 620 ON STATIGN 
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SURFACE SALINITY AND TEMPERATURE OBSERVATIONS 


CRUISE REFERENCE NUMBER 76~ 
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ABSTRACT 


Vertical excursions of instruments suspended beneath a subsurface buoy, 
resulting from hydrodynamic drag forces on the cable, cause problems in 
oceanographic data acquisition and analysis. In an effort to alleviate the 
problem, a plastic cable fairing was evaluated in a field experiment. An 
empirical value for the fairing drag coefficient was obtained using a 
mathematical simulation of the mooring. Then, again using the mooring model, 
comparisons were made of instrument excursions resulting from the use of 
unfaired cable, faired cable and additional system buoyancy. 
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Introduction 


Hydrodynamic drag forces on cables used for mooring an array of 
current meters or other instruments, where the array is supported by a 
Subsurface buoy, result in excursions of the instruments over some interval 
of depth as the current varies in strength. These excursions are undesirable 
inasmuch as they prevent measurements from being obtained at the desired 
vertical positions at all times and they make the task of correlating 
measurements from different stations more difficult. Figure 1 shows such a 
case, where sequential samples from current meter sensors show the vertical 
position vs the current speed at that position over a period of about 2 hours. 
In this particular case, the current meter was susnended 60 m below a 
spherical buoy whose rest position was at a depth of 40 m. Total length of 
%-inch diameter wire rope mooring line was 275 m. Note that the maximum 
vertical displacement shown is equivalent to about 120 m in depth as the 
current speed approached 5 knots. The numbers used to denote the data points 
in Figure 1 refer to direction octants. 


One means of reducina the effect of the cable draq is to increase 
tie buoyancy of the system. Another possibility for improving the situation 
is to reduce the actual drag by fairing the cable to produce a more stream- 
lined profile. An extruded plastic fairing for %-inch diameter cable is made 
by Fathom Oceanology Limited of Port Credit, Ontario, marketed under the 
trade name of "Rigstream" fairing. It was decided to try a field experiment 
to evaluate this product in actual use. A contract was let to Dobrocky 
Seatech Ltd. of Victoria, B.C. to deploy and retrieve a test mooring utilizing 
three current meters on a faired cable and, subsequently, to provide the 
Institute of Ocean Sciences, Patricia Bay with the resulting data tapes for 
analysis, together with a report describing the mooring activities (Dobrocky 
Seatech Ltd., 1977). The drag coefficient for the fairing was unknown, but 
was to be estimated by using the data in conjunction with a numerical computer 
model (Bell, 1977). 


Field Program 


The site chosen for the test mooring was in Haro Strait at the 
position 48°35.6'N, 123°15.5'W. Prior experience suggested that the velocity 
profile at this location was quite uniform in both amplitude and direction, 
especially during periods of large tidal range. The condition of uniformity 
over depth facilitates subsequent analysis. The water depth at the site is 
approximately 720 ft. 


The mooring was deployed on October 22, 1976 and remained in 
position until November 1, when it was retrieved. Because of some handling 
difficulties during installation of the system, the exact amount of fairing 
applied to the cable was not recorded. The fairing segments were each 2 feet 
long and it was known that a few segments were left off so that the personnel 
involved could grip the cable more easily to restrain its motion during 
deployment. It was also noticed, as the cable moved overboard, that a few 
segment of fairing were improperly installed. It is suspected that these fell 
off the cable shortly thereafter. Upon retrieval, the remaining length of 
fairing was obtained and corresponded to approximately 93% of the total wire 
rope length. 
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SCATTER DIAGRAM - PRESSURE VS. VELOCITY 


Figure 1. Vertical excursion of a current meter. System buoyancy - 800 1b, 
Tength of line between current meter and anchor - 215 m, sample 
interval - 10 min, tidal range - 2.5 m. 


The current meters used in the moored array were Aanderaa Model 
RCM-4's, equipped with pressure sensors. The pressure sensor accuracy is 
quoted as being better than 1% of the range. The sensor in the uppermost 
meter of the array had a range of 0-200 psi, with an accuracy equivalent to 
+4.5 ft of saline water. The sensors in the other two instruments had a 
range of 0-500 psi, with a correspondingly reduced accuracy. A further 
consideration is the resolution provided by one data bit in the data logger, 
which is equivalent to 0.5 ft of water. Thus, an uncertainty of one bit in 
each pressure reading could result in an error of +1 ft in determining depth 
differences. A value for the hysteresis error of the pressure sensors is not 
given but, for small changes in depth, it could result in a substantial 
percentage error in depth differences. The accuracy of directions obtained 
from the magnetic compass is given as +5° by the manufacturer. The direction 
and pressure data are obtained as instantaneous values, whereas the speed is 
integrated over the sample period. The sampling interval used in the present 
instance was 2 minutes. 


A spherical steel buoy with a diameter of approximately 3 ft and 
a weight in air of 322 1b was used as a subsurface float for the array. It 
was planned for installation at a height above bottom of about 520 ft, using 
1-inch diameter 6x19 steel cable to connect it to an anchor clump of 3 railway 
wheels and an 80 1b Danforth anchor. The three current meters were to be 
positioned at 100 ft, 300 ft and 500 ft above bottom in quiescent water, i.e. 
a spacing of 200 ft between instruments. The actual separations, as obtained 
by measurement of the components involved and allowing for some cable stretch, 
were slightly different and are given in Table I. The current meters are 
actually installed in the cable by means of a spindle and gimbal arrangement, 
with the cable fastening to the ends of the spindle. The pressure sensor 
location is about 1 ft below the upper end of the spindle, whose total length 
is about 2.4 ft. This is reflected in the entries of the table. 


Sensor Data 


The data records obtained from the test mooring were carefully 
examined to obtain a good sample set for subsequent processing. Unfortunately, 
the instrument in the middle position had an intermittent malfunction which 
narrowed the choice of samples considerably. The chosen data set is presented 
in Table II. There is good correspondence of the velocities. The directions 
are slightly different at the bottom sensor, but not enough to be of much 
consequence. For subsequent use in the mathematical mooring model, the mean 
of 10 samples is used. This is an aid in smoothing out irregularities in the 
data arising from the previously mentioned difference in sampling technique 
between speed and pressure or direction, or from nonuniformities in the 
velocity profile between measurement stations. It should also be pointed out 
that, while oceanographic variables are normally reported in the metric 
system as in Table II, they will here be converted to engineering units (ft- 
lb-sec) for use in the mooring model. These units are used in the model 
because of the present utility of the force unit in connection with mooring 
systems (as anyone who tries to purchase a cable with a breaking Strength of 
X newtons will soon find out). Table III presents some pertinent data thus 
converted. 


TABLE I. Mooring Component Dimensions 


Component 


Anchor clump 
Unstretched cable 
Cable extension 
Swivel & shackles 
Lower spindle end (LSE) 
to bottom pressure sensor 


Bottom pressure sensor 
to upper spindle end (USE) 
Unstretched cable 
Cable extension 
LSE to middle pressure sensor 


Middle pressure sensor to USE 
Swivel & shackles 

Unstretched cable 

Cable extension 

LSE to top pressure sensor 


Top pressure sensor to USE 
Swivel & shackles 
Cable 


Length (ft) 
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TABLE III. Sensor Data Converted To Ft-Lb-Sec System 


Instrument No. 


ralieh!) 2326 2194 


Mean Depth (ft)* 240.9 442.8 624.4 

Rest Position (ft)* 2leon 416.2 616.4 

Change in depth (ft) 28.5 26.6 8.0 
Mean speed = 4.8 ft/sec 


Mean direction = 151°T 


* Depth is obtained from the pressure reading and includes the effect 
due to the water column being saline (specific gravity = 1.024). 


The Mooring Model 


Mathematical simulation of the mooring system is accomplished by 
writing the appropriate static force balance equations for the buoy, the cable 
and the instruments, and solving them simultaneously using a numerical 
integration procedure. Boundary conditions are established by the buoy at 
the upper end of the cable and the integration proceeds from there, in 
discrete steps, to the anchor end of the cable. Here, the anchor and the 
bottom must coincide (within a specified limit) or an iteration is necessary. 
Because of the step-wise integration, the cable consists of an integral 
number of segment lengths and instruments can only be inserted between 
segments. As a result, the model dimensions must be altered slightly from 
those of the real system it is intended to simulate. The effect of this is 
Shown in Table IV for the present case. The next requirement is one of 


setting out the remaining model input parameters with reasonable accuracy, as 
presented in Table V. 


A value for the drag coefficient of the fairing is assumed and the 
model predictions obtained for the change in depth of the sensors under the 
influence of the given current. When this is done for two or three different 
values of drag coefficient, an interpolation can be made for that coefficient 
which corresponds to the measured depth excursions. Mention should be made 
of the fact that only the normal drag force is considered here, consisting 
principally of form drag, but including a small additive skin friction term. 
Longitudinal drag is neglected as an unnecessary complication (although 
provision for calculating it is included in the model). Thus, if the cable is 
at an angle to the flow, the drag is taken basically as the product of the 
square of the sine of the angle and the drag which would occur normal to the 
flow. So long as the angles involved are reasonably close to 90°, the error 


due to the omission of frictional drag in the axial direction is small. 


TABLE IV. Mooring Component Rest Position Distances (ft above bottom) 


Distance 
Component Planned Measured Pressure* Model 
Bottom sensor 100 101.0 Lod 100 
Middle sensor 300 302.4 307.3 300 
Top sensor 500 Uae Uae | 506 
Buoy 520 528.4 528.4 528 


Model cable segment length = 2 ft 
No. of segments = 264 | 


* Derived from pressure sensor data assuming a water depth of 717.5 ft 
so that the top sensor position coincides with the measured value. 


Fairing Drag Coefficient 


Two model simulations were made, using Cy values of 0.45 and 0.55. 
F 


These resulted in predicted depth excursions which bracketed the actual 
excursions of the upper. two current meters. A linear interpolation was then 
made, indicating coincidence of the measured and predicted excursions at a 
value of Cy = 0.50 for the upper instrument, Cy = 0.52 for the middle 
f ‘e 
instrument, and Cy = 0.28 for the lower one. This latter value is suspect in 
F 


light of the agreement in results for the top two meters. Also the bottom 
sensor underwent a relatively small excursion with a possibility of a relative- 
ly large hysteresis error. Much more confidence is placed in the result 
obtained from the top sensor, since it is the most accurate of the three used 
and it experienced the largest excursion in depth. Thus, the value used 
henceforth for the drag coefficient of the fairing (based on frontal area) is: 


Cc. = 0.50 
De 


for sub-critical flow. The uncertainty in this value due to the data logger 
accuracy of +1 bit is only +0.02. More uncertainty resides in the lack of 
complete knowledge of the velocity profile, but a number cannot readily be 
assigned to this. 


TABLE V - Model Input Parameters 


1. Velocity Profile: 
Uniform, 4.8 ft/sec 


2. Buoy: 


Radius = 1.542 ft 

Reynolds No. = 9.9x10° (super-critical ) 

CDp = One 

Buoyancy = 984 lb (for specific weight of water = 64 1b/ft3) 

Weight = 322 1b 

Weight, including the weight in water of swivels and shackles = 356 Ib. 


3. -Cabhe: 
(a) Wire Rope: 


Diameter = 0.021 ft 

Reynolds No. = 6.7x102 (sub-critical) 

Cp = 1.4 

Weight = 10.0 1b/100 ft 

Buoyancy = 2.2 1b/100 ft 

Unstretched length = 515.5 ft (excluding swivels and current meter 
spindles) 


aa 
oO 
al 


Fathom Rigstream Fairing: 


Width (frontal aspect) = 0.030 ft 
Chord *=* 0OS0Z252Ft 

Reynolds No. = 4.0x10* (based on chord) 
Cpe = To be determined 


Weight = 10.2 Ibs/100 ft 
Buoyancy = 7.5 1bs/100 ft 
Length = 478 ft 


(c) Fairing - Wire Rope Combination: 


Width = 0.030 ft 
Weight = 20.2 1bs/100 ft 
Buoyancy = 9.7 |lbs/100 ft 


4. Instruments (Aanderaa current meters): 


Diameter = 0.42 ft 

Reynolds number = 1.3x10° (sub-critical) 

CD, foal 

Weight = 54 Ib 

Buoyancy = 17 1b 

Frontal area = 0.5 ft2 

Instrument locations = 24, 228, 428 ft from buoy 
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At first glance, the estimated drag coefficient value seems high 
for a reasonably streamlined shape. However, it should be emphasized that 
the flow regime is sub-critical, i.e. the boundary layer has not yet become 
turbulent with a consequent rearward movement of the flow separation points. 
(Also, we have usually had much more exposure to examples of super-critical 
drag coefficients for streamlined bodies.) Hoerner (1958, p. 13-19) shows a 
fairing having a cross-section almost identical to that of the Rigstream 
fairing, except that the trailing edge is less blunt. He gives a sub- 
critical drag coefficient of 0.73 for that fairing, and a super-critical 
coefficient of 0.16. Almost always, the cable of any normal mooring system 
will experience only sub-critical flow; towed systems may be subjected to 
Super-critical flow. 


As a check on the validity of the estimate for C the derived 


D b) 
F 
value is used for a prediction of instrument excursion corresponding to 
another point in the empirical data set. At this point, the data is not 
constant over more than 2 or 3 samples, but it is smoothly varying so we can, 
with some confidence, use specific records for a spot check, dispensing with 
averages as used previously. For this case, the converted data is given in 
Table VI. The velocity profile is obviously not uniform but does appear to 
be almost linear. Therefore, a least-mean-squares fit was used to obtain 
coefficients for a profile of the form V = a + bz, where V is the velocity 
and z is the depth. The predicted instrument excursions are also shown in 
Table VI. It is seen that the uppermost two instruments were actually 
displaced downwards more than the distance predicted by the model, with the 
disagreement being of the order of 10%. The indication is either that the 
estimated drag coefficient for the fairing is too small or, more likely, 
that the velocity between instrument positions was higher than assumed. 


TABLE VI. Data Pertaining to Check on Cy Estimate 


F 
Instrument No. 

2199 2326 2194 
Current Speed (ft/sec) 6.2 5.4 Ay? 
Current Direction (°T) 150 157 144 
Recorded Instrument Depth (ft) 269.5 463.6 631.2 
Rest Position (ft) eer 416.2 616.4 
Actual change in depth (ft) 57.1 47.4 14.8 


Predicted change in depth (ft) Biles 43.2 19.6 
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Effect of Faired Cable 


Now that an estimate of fairing drag has been obtained, the 
reduction effected in the magnitude of instrument excursions due to use of 
the fairing can be determined by having the mooring model Simulate both 
faired and unfaired cables at various current speeds. Since the test mooring 
was typical of a normal data-gathering array, it will continue to be used as 
an example. 


Table VII presents the results of several such simulations. An 
examination of this data indicates that a substantial reduction in vertical 
excursion depth is obtained by using the fairing, with the largest 
percentage reduction occurring at the lower velocities. 


TABLE VII. Predicted Instrument Excursions at Several Current Speeds Using 
Faired and Unfaired Cable (All values are in feet). 


Instrument 
Position 


V = 8 ft/sec 
Faired Unfaired 


Faired Unfaired Faired Unfaired 


Top 


Middle 154 


Bottom 


Effect of Increased Buoyancy 


It is instructive to examine the reduction in excursion depth 
which results from an increase in the buoyant force at the upper end of the 
cable. The most practical way of adding buoyancy is to do it in modules 
corresponding to the available buoys, i.e. install additional buoys in the 
system. In practice, this may require the use of a heavier anchor as well. 
If carried to extremes, a stronger (and thicker) cable would be necessary, 
leading to some additional drag, and this might well invoke the law of 
diminishing returns. For the present example, the provision of two identical 
buoys in the system does not require any increase in cable strength. The 
predicted excursions in such a case are given in Table VIII. The results 
there can be compared with the two situations presented in Table VII, where 
it will be seen that a much greater reduction in excursion depth is obtained 
through the addition of an extra buoy as compared to the addition of cable 
fairing. This is more clearly seen, perhaps, in Figure 2, where the data 
pertaining to the top instrument of the moored array is plotted for each of 
the three situations discussed above. 
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TABLE VIII. Predicted Instrument Excursions at Several Current Speeds Using 
Two Buoys and Unfaired Cable (All values are in feet). 


Instrument 3 x - 

Boetetan V = 4 ft/sec V = 6 ft/sec V = 8 ft/sec 
Top 8 Sie) 85 
Middle 7 31 rie 
Bottom 3 14 Oc 


Conclusion 


The preceding considerations show that a cable fairing can be used 
to reduce cable drag and, thereby, the vertical excursions undergone by 
instruments on the cable when the current speed increases. However, a more 
effective way of reducing the excursions is to double the buoyant force at 
the upper end of the cable. In all likelihood, it will be easier to cope 
with an extra buoy during the installation of an array than it will be to 
handle a cable covered with plastic fairing. Some consideration must be 
given, though, to the increased anchor weight required by the use of an 
additional buoy, and to the possibility that a heavier cable might also be 
required if the extra buoyancy increases the cable load beyond the normal 
working range. 


Another point, not yet discussed, concerns the "strumming" 
behaviour of a cable. Strumming, or vibration, results when vortices are shed 
from an object at a frequency close to the natural frequency of the object. 
This is apparently a common occurrence with wire rope and may result ina 
substantial increase in drag. Vortex shedding is inhibited by the use of a 
splitter plate or fairing, so this may well be a consideration in favour of 
the plastic fairing. There is an easier-to-use alternative here as well, in 
the form of the so-called "haired" fairing. This consists of a multitude of 
flexible plastic threads or ribbons fastened permanently to a cable in such a 
fashion as to trail downstream from the cable, acting as splitter plates to 
inhibit vortex shedding. 
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This is a manuscript which has received only limited 
circulation. On citing this report in a bibliography, 
the title should be followed by the words "UNPUBLISHED 
MANUSCRIPT" which is in accordance with accepted bib- 


liographic custom, 


ABSTRACT 


Anomalous oxygen and salinity data from the deep water at Ocean 
Station P (50°N, 145°W) from 1960 to 1977 that were either intentionally 
omitted from publication or were published anyway, fall into two main 
categories - samples exhibiting both low oxygen and low salinity (determined 
by comparison to the average values at that depth or by reversal in the 
sign of the gradient); and samples with low oxygen and 'normal' salinity. 
The low oxygen-low salinity anomalies are attributed to leaking or mistripped 
sampling bottles; the low oxygen-normal salinity anomalies appear to be 
caused by technical problems although no satisfactory explanation is found 
to support this. 


A decline in quality control may be partially responsible for the 
increase in the number of unexplained anomalous oxygen values in recent 
years. 


I suggest that only very obviously spurious data be deleted 
from publication and that all other suspected bad data be published, but 
marked in some way as suspect. 
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INTRODUCTION 


Oceanographic data have been collected at Ocean Station P in the 
North Pacific (50°N, 145°W, depth 4200 m) on a regular basis since 1956. 
Observations at depths greater than 2000 m began in 1960; prior to 1960, 
lack of a suitable winch limited observations to above 2000 m. The data 
were initially published by the Fisheries Research Board of Canada for the 
period 1965 to 1969, [except for 1967 and 1968 when they were published by 
the Canadian Oceanographic Data Centre (now called the Marine Environmental 
Data Service) ], first as the Manuscript Report Series, then as the Data 
Record Series; and were published by the Institute of Ocean Sciences, 
Patricia Bay! inits Pacific Marine Science Report series. Examination 
of the data published for the period of deep water casts, 1976, showed 
occasional deletion of dissolved oxygen content (hereafter termed oxygen) 
and salinity values - presumably because the data was suspect for some 
reason. The increase in the incidence of deleted data by almost 100% in 
1973 drew attention to the increase in the number of irregular values of 
oxygen and salinity occurring, especially in the deep water. Because 
no clear set of criteria had been established for eliminating data before 
publication, we considered it useful to review the original data logs in 
order to re-examine the causes for the frequent data omissions, and at 
the same time, to study any other anomalies present in the data. The 
object of the present study is to determine the reason for the irregular 
oxygen and salinity data in the deep water at Ocean Station P and to 
Suggest a set of criteria for omitting data in the future. 


l 
Formerly known as the Marine Sciences Branch and then the Marine Sciences 


Directorate. 


THE DEEP WATERS OF THE NORTH PACIFIC OCEAN 


The main features of the water column at Ocean Station P can be 
observed in the curves in Figure 1, showing temperature-depth, salinity- 
depth, oxygen-depth and sigma-t-depth for October 9, 1971 - cruise 6, cast 32 
(P-71-6-32). These curves are typical for the time of year. Tabata (1961, 
1965) described the water column as three layered - the upper, seasonal zone 
(0-100 m); the principal halocline (100-200 m); and a lower zone beneath 
the halocline. These layers can be identified in Figure 1. Beneath the 
upper 30 m of isothermal water at 10.5°C, a thermocline is situated between 
30 and 100 m in the seasonal zone. Throughout the seasonal zone, the oxygen 
and salinity are uniform at about 7.5 ml/1 and 32.7%o0. In the principal 
halocline, the temperature is fairly constant although oxygen decreases 
rapidly with depth. Below the halocline, the temperature decreases gradually 
with depth to a minimum of 1.5°C at 4000 m. At 600-900 m, the oxygen reaches 
a minimum of less than 1 m1/1 due to biological decay, after which it 
increases Steadily to about 3.3 ml/1 at 4000 m. Salinity increases gradually 
from the halocline to approximately 34.690%o at 4000 m. The present study 
is concerned only with the characteristics of the deep and bottom water below 
3000 m. 


The water below 2500 m at Ocean Station P is part of the Pacific 
Deep and Bottom Water called Oceanic Common Water by Montgomery (1958). Its 
flow northward from the Antarctic Circumpolar current was inferred by Knauss 
(1962) from the distribution of water properties and further studied by 
Reed (1969), but the circulation has been difficult to determine because of 
the small variation in horizontal water properties, and is probably far more 
complicated at smaller spatial scales. Fofonoff and Tabata (1966), for 
example, detected a relatively narrow "tongue" of high salinity water 
extending from the north in the deep water at Ocean Station P. 


Reed (1969) studied average values of water properties in the North 
Pacific at 3000 m, 4000 m, and 5000 m based on United States Coast and 
Geodetic Survey Cruises from 1961-1966, between 20°N and 60°N, and 150°W and 
170°W. At 4000 m, the range of salinity over the whole region was 34.675%0 
to 34.700%0 and the range of oxygen was 3.30 ml/1 to 3.60 ml/1]. At the 
reoccupied sites, Reed found that 95% of the salinity values deviated less 
than 0.008%o from the mean Salinity over 5 years; and that 95% of the oxygen 
values deviated less than 0.12 ml/1 from the mean oxygen value. During this 
particular 5 year period then, the variation in deep water properties. over the 
entire region was small. At Ocean Station P itself, Rupp (1969) found the 
average properties at 4200 m near the bottom from 1961-1965 to be 1.50°C, 
34.69/00, and 3.30 ml/1 and their respective standard deviations to be 0.02°C, 
0.012 %0 and 0.10 ml/1. 


Figures 2a and 2b show six T-S and T-05 curves that are fairly 
typical of the water from 3000 to 4200 m at Station P. In general, as the 
temperature decreases with depth, both oxygen and salinity increase. 
Occasionally, there are small decreases in salinity (as in P-62-3-13 by 
0.005%0o, Figure 2a), or in oxygen, but this is not unusual considering the 
small gradients at 4200 m and the limits of the errors (+0.003%o and 
+0.03 ml/1) defined by the methods of analysis described in the next section. 


There are also increases in temperature of up to 0.03°C with depth near the 
bottom as seen in Figures 2a and 2b - eg. 0.02°C in P-60-2-19. 


ANALYSIS OF THE DATA 


1. Original Measurements and Computations 


All temperatures given in the present study are observed 
temperatures. Since 1969, they were recorded using a pair of Richter & 
Wiese reversing thermometers on each sampling bottle. The thermometers 
yield an accuracy of +0. 02°C in routine use, and the accuracy based on 
readings from two thermometers should be greater. 


Salinities were measured using an inductive salinometer--ashore 
from 1960-1962, and in the ships' laboratory from 1962 to date. The two 
models presently in use are an Auto-Lab Model 601 MK III and a Plessey 
Model 6220. Both salinometers are calibrated with Copenhagen water approx- 
imately every 20 samples. The accuracy using duplicate determinations is 
estimated to be +0.003%o in a shore laboratory and +0.004%o in a ship's 
laboratory at sea (Strickland, 1958). 


The dissolved oxygen content was measured using the Winkler or 
modified Winkler method with an estimated error of +0.03 ml1/1 (Strickland 
and Parsons, 1968). 


Two types of sampling bottles were used to collect the water 
samples - Nansen bottles from 1960 to 1972 and Niskin bottles from 1972 to 
1976. The practice in drawing both salinity and oxygen samples is to first 
rinse the bottles and caps twice; then a standard method is used to draw 
each sample although individual techniques may vary. 


2. Omission of Data from Publication 


First, all omitted data points were checked against the original 
data logs and when found, recorded; then an attempt was made to determine the 
approximate guideline used to eliminate this data. The oxygen values deleted 
fell into two categories - either the absolute value of the oxygen was too 
high (or low) compared to the average value for that depth (eg. greater than 
3.6 ml/1 or less than 2.6 ml/1 at 4200 m) or there was a reversal in the 
vertical gradient of oxygen, causing an apparent decrease in oxygen with 
depth. In most cases, the oxygen had to decrease at least 0.1 ml/1 from one 
sampling bottle to the next deeper one for the value to be omitted (aithough 
there were cases of oxygen values being omitted when the difference was as 
low as 0.07 ml1/1). Salinity values were omitted when a decrease in salinity 
of greater than 0.02%o between successive sampling bottles occurred. 


Table I shows the number of oxygen and salinity values omitted each 
year from the deep water data for unexplained irregularities. The total 
number of deep casts taken each year is recorded in Column 1. Column 6 shows 
that from 1967 to 1973, 25 oxygen values were omitted for no apparent reason, 
and Column 7 shows that 26 oxygen values and 35 salinity values were omitted 
for documented technical reasons during the same period. Technical reasons 
included: 

1) bad oxygen reagents - 5 cases, 

2) overtitration of oxygen sample - 3 cases, 


Table I. A Table of Oxygen and Salinity Values Omitted From Publication 
From 1960-1976 


Year No. of No. of No. of No. of No. of No. of 
Irregular | Irregular | Irregular | Irregular| Values Values 
02 Values| S%o 0. Values| S%o Omitted | Omitted 
Published | Values Omitted Values For No For 
Published Omitted | Apparent | Technical 
Reason Reasons 
1960 92 $%o 
1961 - 
1962 - 
1963 - 
1964 ae 
1965 /4 
1966 /| 
1967 = 
1968 /5 
1969 /| 
1970 2/5 
1971 4/3 
1972 8/3 
1973 5/4 
1974 =/4 
1975 1/5 
1976 4/1 
Column #: 1 2 3 4 5 6 7 


) bottom sediment in sampling bottle - 1 case, 

) leaking sampling bottle - 10 cases, 

) drift in salinometer - 15 cases, 

) nonrepeatability of a salinity value - i.e. duplicate samples 
more than 0.01%o apart - 7 cases, 

7) Others - 9. 


Next, the published data was examined for irregular data points that 
had not been omitted. Columns 2 and 3 show the number of irregular oxygen 
and salinity values that had been published. The criteria for selecting 
these values as irregular were as follows: 
) oxygen decreasing with depth by greater than 0.01 ml1/1 between 
sampling bottles, P 
2) salinity decreasing with depth by greater than 0.0120 between 
sampling bottles, 
3) an oxygen value that was at least 0.3 ml/1 higher than the 
average value for that depth. 


Although the data has a published accuracy of +0.04%o for duplicate 
samples of salinity run at sea and 0.03 ml/1 for oxygen, it is felt that 
the above uncertainties (+0.05 ml/1 and +0.006°%0) are more realistic as 
error limits for testing irregular data. For example, the published errors 
do not take into account the variations in drawing the samples before the 
determinations of oxygen and salinity are made and the possible errors 
introduced at that time. The published errors assume the samples are run by 
highly trained technicians which may not always be the case. 


Rupp (1969) showed that in the deep water at Station P from 1960- 
1966 the actual uncertainties encountered were 0.012 %0 and 0.10 ml/1l. The 
Study by Reed (1969) for the same period at several other North Pacific 
Stations showed a deviation from the mean of less than 0.008%o for Salinity 
and 0.12 ml/1 for oxygen in more than 95% of the cases out of 85 samples. 
Corritt and Carpenter (1966) showed that random and systematic errors in 
oxygen data using Winkler titrations may be as great as 0.1 ml1/1 and that for 
Saturated samples and very small values the errors could be much greater. 
Strickland and Parsons (1968, page 21), in their description of the Winkler 
method, state that the precision they predict; 0.03 n ml/1, where n is the 
number of samples run, is "the highest precision considered likely for work 
in a shorebase laboratory under near ideal conditions, using thiosulphate 
Standardized by the mean of at least five titrations" and that "under routine 
conditions at sea, the uncertainty range will be appreciably greater, 
perhaps doubled". 


A recent, unpublished study of the well mixed surface layer at 
Ocean Station P by S. Tabata shows that from 1969 to 1976 the mean standard 
deviations of the water properties in the well] mixed surface layer were as 
follows: 


October-March April-September 
Temperature +0.02°C +0.03°C 
Salinity +0.003%o +0.003%o 


Oxygen +0.04 ml1/1 +0.04 ml/1 


The above statistics are drawn from well over 100 observations. 
They show the mean standard deviation of salinity values to be the same as 
the predicted error of +0.003%o and the mean standard deviation of the 
oxygen values to be slightly larger than the published error of +0.03 ml/1. 
The mean standard deviation in the temperatures is slightly larger in the 
period April to September than the period October to March, probably due to 
the direct heating effect of the sun. 


The number of published, irregular oxygen and salinity values 
recorded in Table I demonstrates the arbitrariness of data omissions. For 
example, in 1975, an extreme case, 17 irregular oxygen values were omitted 
from publication but 10 irregular oxygen values were not omitted. For most 
years (Table I), there are a few irregular oxygen and salinity values not 
omitted, even though they appear to fit the estimated criteria previously 
used. There are also several blocks of data which appear suspect but which 
have been published anyway. The three examples are: 

1) P-76-1: all the temperatures from 3500 m-4200 m are 0.03 to 

0.06°C higher than normal for those depths. This has been 
attributed by C. de Jong, chief technician for Ocean Station P 
cruises (personal communication) to uncalibrated thermometers, 

2) P-61-5: Casts 8 and 13 - low salinities encountered in all the 

deep samples; eg. 34.653%o0 and 34.664%o at 4200 m, 
3) P-76-6: low salinities in all samples from 3000-4200 - 
eg. 34.662%o0 at 4200 m. 


The following is a list of data that could not be found with the 
appropriate original log sheets. Either the data were misfiled, mislaid or 
were never computed. 

1) P-63-1: Cast 12; oxygen values not published and original data 

Sheet not found. 


2) P-65-2: Cast 11; oxygen values not published and original data 
not found. ? 

3) P-70-4: Casts 3 to 9, no deep oxygens calculated. 

4) P-70-5: Cast 46; the computations were not completed. 


Finally, it should be mentioned that there are many anomalous 
Salinity values omitted from the 1976 data reports and replaced with inter- 
polated values but are not marked as interpolated in these reports. Also, 
some of the published salinity values are the 'better' value of two duplicate 
samples run from the same sampling bottle, even though the two values are 
separated by more than 0.01%o0. The usual practice is to run duplicate 
salinity samples from each depth, then to average the two values obtained to 
arrive at the published salinity value. If the difference between the two 
values is greater than 0.01%o then the values should be discarded and an 
interpolated value inserted. 


3. Significance of the Irregular Data 


When large amounts of data are collected over a period of years such 
as at Ocean Station P, there are certain to be spurious data points, even 
with strict quality control. Over the 17 years under study, 50 different 
people have taken observations at Station P, some with no previous deep sea 
experience and little training in the observational and analytical procedures 


used to obtain oceanographic data aboard ship. Therefore, it is not 
unexpected to find a few anomalous oxygen or salinity values in the data. 
After 1967 though, the number of irregular oxygen and salinity values 
increases so markedly that it is apparent these errors are not merely random. 
Table II lists the total number of irregular oxygen and salinity values each 
year as defined by the criteria on page 6. It also lists the number of 
anomalous oxygen and salinity values per deep cast each year. Note the 
increase through the years. Up until 1968, except for 1963, there were fewer 
than 1 anomalous oxygen or salinity value for every 5 casts. In 1963, 3 

of the irregular oxygen values were high values associated with one particular 
cast. From 1968 to 1974, the average number of irregular oxygens increased 
to 1 every 3 casts; the average number of irregular salinities also increased 
to 1 every 3 casts. By 1975-76, there was more than one irregular oxygen 
value per cast and about one irregular salinity value per cast. 


The most common recurring anomaly is combined low oxygen-low 
salinity. Prior to 1968, this event never occurred. Another, less common 
anomaly is a low oxygen value accompanied by a 'normal'! salinity value. 
Except for a few cases, the temperature readings encountered are not anomalous. 
The following is a summary of low oxygen-low salinity; and low oxygen-'normal ' 
Salinity occurrences from 1968 to 1976. 


1968 P-68-1 - every deep cast had low oxygen, low salinity in the 4200 m 


bottle. 
1969 - irregularities sporadic 
1970 - irregularities sporadic 


197] P-71-1 - 3 low oxygen-low salinity out of 5 casts (one documented as 


leaking). 


P-71-3 - every deep cast exhibited low oxygen-low salinity in the 4000 
m (next to bottom) bottle although two of these cases were within 
limit of error. 


1972 P-72-3 - low oxygen in the bottom bottle for two out of three casts. 


P-72-/7 - three of five deep casts had low oxygen-low salinity in the 
bottom bottle, 2 out of five deep casts had low oxygen-low salinity 
in the 3000 m bottle. 


- every deep cast (4) had low salinity-low oxygen in the bottom 


Starting in 1973, the 4200 m bottle was replaced with 2 bottles 
10 m apart. This was done specifically to investigate the already increasing 
irregularities observed in the deep water. 


1973 P-73-9 - in every one of 5 casts, at least one of the bottom bottles 
had Tow oxygen-low salinity. 


1974 P-74-7 - low oxygen-low salinity for casts 717, 721, 730 - twice in a 
3000 m bottle. 


1 positive gradient with increasing depth, and within 0.012%0o of the 


average salinity for that depth 
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P-74-9 - casts 923, 911, 917, 930 - low oxygen-low salinity in top 
bottle at 4200 m. 


1975 P-75-1 - casts 109, 120, 127, 133, 134 - low oxygen-low salinity in 
one or both bottom bottles in each cast. 


P-75-3 - casts 309, 320, 323 of 6 casts - low health low salinity in 
bottom bottle. 


P-75-4 - no irregularities. 


P-75-5 - no consistent irregularities. 


P-75-7 - low oxygen-low salinity in four of six casts 
low oxygen only in two of six casts 


P-75-9 - low oxygen-low salinity in 3 out of 6 casts 


1976 - very low oxygen values (1.4-2.4 ml/1) and low salinity 
ee a 34.5%o0) at 4000 m for all four casts. 


P-76-2 - relatively low oxygen values in six out of six casts in 
bottom 4200 m bottle (2.8-32. ml1/1). 


P-76-4 - 5 cases of very low oxygen-low salinity in bottom 4200 m 
bottle; 4 other cases of low oxygen-low salinity occurred. 


P-76-6 - 6 deep casts; 2 cases of low oxygen-low salinity in both 
bottom bottles; 4 cases of low oxygen-low salinity in only one bottom 
bottle. 


This summary shows a marked increase in 1975 in the number of 
cruises where irregular values occur consistently. It also shows, that 
looking at the whole period from 1968-1976, the occurrence of cruises 
recording recurring low salinity-low oxygen is really quite sporadic and 
does not seem to show any kind of a pattern. Most of the odd values are found 
in one of the bottom two bottles (at 4200 m). From 1968 to 1976; of the 70 
cases of low salinity-low oxygen, 51 or 70% were found in one of the bottom 
bottles. In P-72-3 and in P-76-2, the anomalies were low oxygen accompanied 
by "normal" salinity. 


The family of T-S and T-0, curves in Figures 3a and b and 5a and b 
illustrates the two anomalous situations encountered. Figures 3a and b show 
consistently low dissolved oxygen in the lower 4200 meter bottle during cruise 
P-76-2. There is no consistent irregularity in the accompanying salinity 
although in several casts (particularly noticeable in cast 2-5) the salinity 
decreases slightly in the bottom bottle. (The T-S curve and the T-05 curve 
in Figures 4 and 4b respectively are drawn from the mean temperature, 
Salinity and dissolved oxygen values for each bottle from six casts and make 
the general trends of the curves clearer.) The mean T-05 curve of Figure 4 
exhibits the same general pattern as the 'typical' T-05 curves in Figure 2b 
except that by oxygen in the bottom bottle decreases by approximately 0.3 
ml/1 from the bottle 10 m above. The corresponding T-S curve in Figure 4a 
exhibits relatively low salinity in the bottom bottle compared to Figure 2a - 
34 .675°%0 compared to 34.680%o but shows no appreciable reversal in the 
salinity gradient. One unusual factor of the large reversal in vertical 
oxygen gradient is that it occurs only between the bottom two bottles which 
are 10 m apart; even though the range of depth of the bottom bottles is 3925 
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to 4183 m - a range of over 250 m (the published error in depth is +5 m). 


Figures 5a and 5b, P-71-3, illustrate the situation of low 
salinity-low oxygen. Again, the mean curves produced in Figures 6a and b 
show the trends more clearly. The unique feature is the decrease in 
dissolved oxygen (0.2 m1/1), temperature (0.02C°), and salinity (0.09%o) in 
the second to bottom bottle at 4100 m (here, the last three bottles are 
each separated by 100 m). During this cast, the second to bottom bottle 
ranged over a depth from 3895 to 4090 m; the bottom bottle actually sampled 
within this depth range during cruise P-71-3 but recorded "normal" salinity 
and oxygen. The decrease in salinity at 4000 m makes the water column 
unstable at this depth. Sigma-t decreases as much as 0.015 over 100 m, 
eg. Wr s7720to 27 2657". 


There are two possible reasons for the recurring anomalies - 
either they are caused by a real phenomenon or by a recurring technical 
problem. One technical problem that is occasionally referred to in the log 
sheets, and could be responsible for the low oxygen-low salinity situations, 
is a leaking sampling bottle. Between 1970 and 1974 there were 10 recorded 
cases of leaking bottles; none were recorded before 1970 and none in 1975 
and 1976. C. de Jong (private communication), states that Niskin sampling 
bottles are susceptible to leakage if the top and bottom lids are not seated 
properly upon snapping shut, especially if this occurs with an open air vent. 
Sometimes, if Niskin bottles are left cocked for extended periods, a nick 
can form in the O-ring which could contribute to leakage. de Jong also 
states that a leaking bottle might not be noticed by an observer, because 
even perfectly functioning sampling bottles drip a lot of sea water under 
normal conditions. 


The characteristic properties of a leaking sampling bottle seem to 
be low salinity and low dissolved oxygen in varying degrees (although 
recorded cases don't always exhibit low oxygen). Table II lists the number 
of cases of low salinity-low oxygen occurring each year. The number increases 
quite considerably in 1972, the year the sampling bottles were changed from 
Nansen to Niskin. It is possible that Niskin bottles are more susceptible 
to leakage. If it is assumed that the low salinity-low oxygen cases are 
explained by this technical failure, then the remaining number of unexplained 
anomalous oxygen values per case each year are reduced to the ratios in 
Column 6, Table II. Except for 1963 (previously discussed), and 1975 and 
1976, most of the ratios for the other years are low enough that the remaining 
anomalies are probably random. In 1975, half (6) of the remaining irregular 
oxygen values after the ‘leaking bottle’ cases are removed, are high values 
scattered throughout the cruises. The rest are low oxygen values, 3 of 
which are found in the bottom bottle of 3 casts during P-75-7 and P-/5-9. 

In 1976, five of the remaining 11 irregular oxygens occur in the lower 4200 m 
bottle during P-76-2. The oxygen values range from 2./ mi/ 12 tew3.dimi/13; the 
oxygens decrease 0.2 to 0.5 ml/1 from the bottle 10 m above. Thorough 
examination of the original data reveals no computation errors to account for 
this consistent decrease in oxygen at 4200 m. As shown in Column 8, the 
number of anomalous salinities per cast is also much lower after those low 
salinities occurring with low oxygens are removed. Even in 1975 and 1976, 
the number of irregularities per cast is not excessive. 
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4. Paired Bottles 


Starting in 1973, another sample bottle was placed 10 m above the 
bottom bottle at the 4200 m depth for every deep cast. The yearly statistics 
for these paired bottles are shown in Table III. The number of cases of 
irregular values is far greater in 1975 and 1976. 


i Columns 2 and 5, Table III, show the mean oxygen (05) and salinity 
(S%o) values for both upper and lower bottles at 4200 m for each year. If 
the average of the two 0, values is calculated for each year, and the yearly 
averages compared for the four year period, a yearly decrease in the average 
0, is observed: from 3.28 ml1/1 in 1973 to 3.12 ml/1 in 1976. This compares 
with 3.3 ml/1 observed by Rupp (1969) as the average deep water oxygen from 
1961-1965. Over this same period, the average temperature rose from 1.50°C 
(Rupp, 1969) to about 1.52°C (average temperature for 1973-1976). Another 
noticeable change is the_increase in the difference between the upper bottle 
0> and the lower bottle 0, from 0.05 m1/1 in 1973 to 0.25 ml/1 in 1976. This 
can also be seen in Figures 7a and 7b, time plots of the individual differ- 
ences in salinity (a) and oxygen (b) between the upper and lower bottles from 
1973 to 1976. The dashed lines show the published error in salinity and 
oxygen values (0.006%o and 0.06 ml/1) that are increased slightly to include 
the error in drawing the samples (to 0.012 and 0.1 ml/1). A positive 
difference indicates a higher salinity or oxygen in the bottom bottle and 
vice versa. The most outstanding feature is the increased number of "A0," 
values outside the dashed error lines in 1975 and 1976 compared to 1973 and 
1974. The larger "AQ>" are more prominent during certain months such as 
December 1974, January to April 1975 and February to August in 1976. Salinity 
differences, "AS%o", are also large for May, June, July and August of 1976. 
Another period of high "AS%o" is August 1974 to January 1975. 


The yearly standard deviation of the means in Table III includes 
Systematic and random errors as well as real yearly variations. The standard 
deviations of the mean oxygen values range from 0.14 to 0.35 ml/1. Reed 
(1969) found that the range of dissolved oxygen in the deep water of the North 
Pacific was 0.40 ml/1, so in comparison, considering this horizontal range of 
properties, the standard deviations are not entirely unreasonable, although 
0.35 ml/1 seems a high yearly standard deviation for 4200 m, especially 
considering that the bottle only 10 m above had a standard deviation of only 
0.13 ml/1 for the same period. The standard deviations of the mean salinities 
range from 0.005%o to 0.06%. Again, the large differences between the 
Standard deviations of the upper and lower bottles seem anomalous. The 
range of salinity found by Reed (1969) in the study mentioned was only 
0.02%o, much smaller than the yearly standard deviations calculated for 1973, 
0.06%o and 0.05%o. 


Column 6 shows the 95% confidence limits of difference between the 
pairs of S%o and 0» values. In general, zero lies between the confidence 
limits indicated except for the salinity in 1974 and the salinity and oxygen 
in 1976. This implies that in 1976, there was a 95% probability that in 
1976 both the oxygen and the salinity were actually higher in the upper 
bottle. 
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There are cases of duplicated anomalies. In each of the years 
1973 and 1974 there was one case of low oxygen in both the upper and lower 
bottle. In both bottles, the oxygen values were within 0.03 ml/1 of each 
other (2.84, 2.85 ml/1 and 2.61, 2.64 ml/1) - well within the published 
error. In 1975, three out of six of the simultaneous cases of low oxygen- 
low salinity were within 0.008%o and 0.05 ml/1 of each other. In one of 
the cases of low salinity in both bottom bottles in 1976, the low values of 
salinity were identical (34.662%o0); and in one of the cases of low oxygen- 
low salinity in both bottom bottles, the salinity values were identical 
although the oxygen values were different by 0.2 ml/1. 
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DISCUSSION 


Many of the anomalous oxygen and salinity values in the deep water 
at Station P appear to be caused by leaking sampling bottles. The 
suspected cases of leaking sampling bottles exhibit the same low oxygen-low 
salinity combination found in the few recorded cases which display water 
properties characteristic of shallower water. The range of salinity, 34.5%o 
to 34.7%0 (excluding salinity values below 34.0%o), and of oxygen, 1.4 ml/1 
to 3.2 ml/1, in these anomalous cases seems too large for them to be 
indicators of a single real phenomenon. The low oxygen-low salinity cases 
each year are usually concentrated in certain cruises rather than scattered 
throughout the year, but the frequency of cases and date of these cruises 
vary from one year to the next - for example, in 1971, there was only one 
cruise in April-May with this anomaly recurring at 4200 m; in 1972, there 
were two cruises - one in October and one in December; in 1973, there was 
only one cruise in December exhibiting low salinity and oxygen at 4200 m; 
in 1974, there were two cruises - one in October and one in November; in 
1975, low oxygen and salinity occurred repeatedly in 3 cruises - at 4200 m 
in January, at 4000 m in October and at 4000 m in December; and finally, in 
1976, up to November, low oxygen and low salinity occurred consistently at 
4200 m in two cruises, one in May-June and one in August-September. Again, 
the cruises with anomalous values seem to occur too sporadically to attribute 
the anomalous salinity and oxygen values to a single phenomenon. 


There are many speculations about the circulation of deep water in 
the North Pacific inferred from the distribution of water properties. Knauss, 
(1962) for example, predicted that the deep water moves up the Pacific Basin 
from the south with a net speed of 0.05 to 0.1 cm/sec and is gradually 
warmed by mixing with water above and from heating from the earth's interior. 
The mixing also gradually decreases its salinity. He believed that the deep 
water moved more readily up the western side and there were indications of 
an eastward flow at the equator. Reed (1969), on the other hand, suggested 
an eastward transport of water along 50°N, south of the Aleutians and north- 
ward transport at 5000 m south of 40°N. Because of the small variations in 
horizontal properties, deep circulation is difficult to infer. Direct current 
measurements in the deep Pacific have proved thus far unreliable because of 
the limit of the current meters used in resolving water speeds down to 1 mm/ 
sec from 0-5 cm/sec (Collins, 1969) and the fluctuating barotropic flows of 
1 cm/sec or more (Barbee, 1965) and deep tidal currents at least this great, 
which make resolution of small (order of 1 mm/sec) net velocities difficult 
(Knauss, 1962). The measurements taken by Collins (1969) using a Geodyne 
Current Meter model 102 with a Savonius rotor on the bottom (4250 m) at 
Ocean Station P recorded mean current speeds of 1.54 cm/sec and 0.47 cm/sec 
on two separate days. The net transport velocities are generally considered 
small. It is impossible to explain the extremely low oxygens and low 
Salinities and large variations in the time scale of days by such slow 
advection, especially considering the small range of properties at 4000 m 
over the whole North Pacific from 1961-1966 found by Reed, i.e. 0.3 ml/1 
and 0.025%o. It is improbable that the instability caused by a large 
decrease in salinity with depth at constant temperature would be present at 
4000 m for any length of time. 
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Three other factors suggest that the anomalous values are due to 
leaking sampling bottles: first, the low oxygen-low salinity samples usually 
come from a particular bottle on each cruise rather than a particular depth, 
which implies that some bottles are more susceptible to leaking than others: 
second, in the 44 cases of low salinity-low oxygen in paired bottles at 
4200 m, from 1973-1976, only 9 occur simultaneously in both bottles and in 
only one case are both the salinity and oxygen of the two depths within 
0.012%o0 and 0.1 ml/1.of each other. If the anomalies are caused by such a 
localized pheonomenon that only one of two bottles 10 m apart can detect it, 
it seems unlikely that the anomalous characteristics would be sampled so 
repeatedly in the same bottle. Third, the temperatures which are measured 
at depth by reversing thermometers show no unusual deviations from the average 
values for the deep water in these low salinity-low oxygen samples. If the 
samples are indicative of a different water mass, it is probable that the 
characteristic temperature would also be different. 


It is possible that some of the low oxygen-low salinity samples 
are due to sampling bottles that do not shut when originally hit by the 
messenger but trip on the way up to the surface due to a sudden jerk on the 
wire. B. de Lange Boom (personal communication), an observer aboard the 
weathership in January, 1976 states that the present messengers are poorly 
suited to the Niskin bottles and occasionally fail to trip a bottle. 
Examination of the cases of the deep water low salinity-low oxygen samples 
from 1971-1976 shows that in 1976, 7 out of 20 of these samples exhibited 
oxygen and salinity corresponding to a specific shallower depth - usually 
the shallower depth is 3000-3500 m; in 1975, 8 out of 17 low oxygen-low 
salinity samples are from a specific shallower depth; in 1974, 2 out of 8 
cases fall into this category; in 1973 - none; in 1972 - 1 out of 10; and 
in 1971 - 3 out of 6. These cases could result from mistripped bottles. 


If leakage or mistripping is the technical difficulty causing the 
low oxygen-low salinity irregularity, one unsolved characteristic of the 
problem remains: why do the anomalies occur so often in the bottom bottle? 
Looking at the water column above 800 m, where salinity is increasing and 
oxygen is decreasing with depth, there are seven cases of low salinity and 
high cxygen recorded in 1975 and 1976, which would be characteristic of a 
leaking sampling bottle at these depths. The leaking bottles are not 
exclusive to 4200 m then, but do occur there much more often than at other 
depths. Also, in the few cases of duplication of irregular values in the 
two 4200 m bottles, it is hard to explain the anomalies by two sampling 
bottles leaking in exactly the same way or mistripping at exactly the same 
depth. 


If the low salinity-low oxygen samples are caused by leaking 
sampling bottles, the only remaining unexplained anomalies are the low oxygen- 
normal salinity samples found in the 4000-4200 m bottles during P-72-3, 
P-75-7, and P-/6-2. The low oxygens occur in anly 2 out of 3 casts in 
P-72-3. In P-75-7 and P-75-9, the 3 low oxygen-normal salinity samples are 
found in the bottom bottle of 3 casts interspersed between 6 other casts 
exhibiting low oxygen and low salinity in the bottom bottle. These anomalies 
in themselves do not appear significant. But in P-76-2, every bottom bottle 
exhibits low oxygen (from 0.2 to 0.5 ml/1 lower than normal) although two of 
these samples also exhibit low salinity as well. These low oxygens appear to 
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be caused by some technical problem: first, because they always occur in 
the same bottle during P-76-2; second because the low oxygen values are not 
reproduced in the sample from 10 m shallower, even though the difference in 
oxygen in the two bottles is as much as 0.5 ml/1 and the range of depth of 
the paired bottles is 250 m. Unfortunately, there is no obvious explanation 
for the anomalous oxygen values. 


Although oxygen is usually considered a conservative property in 
the deep water, it can be changed by biological influences. Sedimentary 
respiration affects the oxygen only a few centimeters above the bottom, 
and would not likely be detected at the depths nor to the degree of the 
recorded oxygen anomalies (Dr. Grant Gardner, Institute of Ocean Sciences, 
Patricia Bay, personal communication). A strong density gradient might trap 
organic material, causing a depletion of oxygen, but there are no such 
gradients associated with these low oxygens. Since 1962, no nutrient analysis 
has been performed on deep water samples so it is not possible to use them 
as indicators of biological activity. Observations of nutrients in 1960 
Show silicate increasing from the surface to the bottom with the highest 
values (app. 162 yg.a./1) from 1000 to 4000 m; the limited number of 
phosphate observations in 1960 show phosphate decreasing from 3 ug.at/1 at 
1500 m to 2 ug.at./1 at 4000 m. 


It should be noted that the anomalous oxygen values are not 
entirely limited to the deep water; for example, cruise P-76-2 recorded 
high (by at least 0.25m1/1) oxygen values in the 800 m bottle in the oxygen 
minimum layer in every cast. 


The low oxygen values in the deep water at Ocean Station P that 
are not accompanied by low salinity in 1975 and 1976 defy explanation as 
either a technical problem or a real phenomenon. If the low oxygen is a 
real phenomenon, the many cases of low oxygen-low salinity caused by leaking 
sampling bottles would tend to mask any low oxygens that were real. Sporadic 
cases of high oxygen-normal salinity, non-duplicated salinity, and low 
salinity-normal oxygen which occur in 1975-76 also make interpretation of the 
low oxygen cases difficult. 


19 


CONCLUSIONS 


Samples taken from the deep water that have both low oxygen and low 
Salinity are attributed to leaking sampling bottles or mistripped bottles. 


The samples taken from the deep water with low oxygen and normal salinity 
during P-76-2 are essentially inexplicable but are probably also caused 
by sampling or analytical technique. 


A firm set of criteria should be set down for eliminating data from 
publication. Probably only recorded technical failures and extremely 
irregular values (especially in the case of samples exhibiting low 
oxygen, low salinity and normal temperature which indicate a leaking 
sampling bottle) should be omitted, with other suspected data being 
published but marked in some way to show they are suspect. The most 
recent reports are beginning to use this procedure, noting both 
interpolated and suspected values. When determining whether values are 
irregular, the published error should be slightly increased to include 
error in drawing the samples, so that two samples are considered 
Significantly different if they are separated by 0.01 to 0.012%o for 
0.08 to 0.1 ml/1, for example. Suspected values should be those a) that 
show reversal in vertical gradient by at least the above amounts; or 

b) oxygen or salinity values that differ by more than 0.3 m1/1 or 0.02%o 
from the average values for that depth (especially if a whole cast 
exhibits anomalous values). All interpolated values should be marked, 
including those inserted to replace omitted original data. 


The increase in the number of unexplained anomalous oxygens, the 

increase in difference between the oxygen values in the two 10 m bottles 
at 4200 m, and the increase in the number of cases of non-duplicated 
salinities in recent years seem to indicate that a higher degree of 
quality control in data collection and reduction should be exercized in 
the future. More care should be taken in recording all technical failures. 
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Figure 1: T-Z, o¢-Z, S%o0-Z, 02-Z plots for Ocean Station P, October 9, 1971. 
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Preface 


This report is intended to serve as a user's manual to G. Godin's 
tidal heights analysis and prediction programs, revised along lines suggested 
by Godin. In addition to describing input and output of these programs, the 
report gives an outline of the methods used; a full presentation of which 


can be found in Godin [1972] and Godin and Taylor [1973]. 


Users who wish to receive updates of these programs and manual 


should send their names, addresses and type of computer used for implementa- 


tion to the author. 
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An USE OF THE TIDAL HEIGHTS ANALYSIS COMPUTER PROGRAM 


wal General Description 


This program analyzes the hourly height tide gauge data for a given 
period of time. Amplitudes and Greenwich phase lags are calculated via a least 
squares fit method coupled with nodal modulation for only those constituents 
that can be resolved over the length of record. Unless specified otherwise, a 
standard data package of 69 constituents will be considered for inclusion in 
the analysis, however, up to 7/7 additional shallow water constituents can be 
requested. If the record length is such that certain important constituents 
are not included directly in the analysis, provision is made for the inference 
of the amplitude and phase of these constituents from others. Gaps within the 


tidal record are permitted. 


1.2 Routines Required 


1) PENG coast cee reads in some of data, controls most of the output and 


calls other routines. 


2) ENPUT act. yd ue reads in the hourly height data for the desired time 


period and checks for errors. 


3) UCON SY 2.2. oo! chooses the constituents to be included in the analysis via 


the Rayleigh criterion. 


4) SCR] HZ. o7.. 83 finds the least squares fit to an equally Spaced time 
series uSing sines and cosines of specified frequencies as 


fitting functions. 


5) MURS vac troraeatoe y reads required information and calculates the nodal 


corrections for all constituents. 


6) INFERS eee reads required information and calculates the amplitude and 
phase of inferred constituents, as well as adjusting the 


amplitude and phase of the constituent used for the inference. 


7) GHLSKYs<-.o204 solves the symmetric positive definite matrix equation 


resulting from a linear least squares fit. 


8) GUA Yi stu. csee returns the consecutive day number from a specific origin for 


any given date in the range 1901 to 1999, and vice versa. 
9) SET TABS. xa <5 sorts results in preparation for ‘J' card output. 


10) DGHS88 2.0. . «: produces ‘J' card output and lists the constituents in species 


by increasing frequency. 


pats: Data Input 


For a computer run of the tidal heights analysis program, two logical 
units are used for data input. Logical unit number 8 contains the tidal 
constituent information while logical unit 4 contains the hourly heights and 


information relating to the type of analysis and output required. A listing of 
the standard constituent information for logical unit 8 and a sample set of 
input for logical unit 4 are given in Appendices 1 and 2 respectively. 


Lodical unit 8 expects four types of data: 


i) One card each for all the possible constituents, KONTAB, to be included 
in the analysis along with their frequencies, FREQ, in cycles/hour and 
the constituent with which they should be compared under the Rayleigh 


criterion, KMPR. The format used is (4X,A5,3X,F13.10,4X,A5). Unless 


chi 


144) 


KONTAB is specifically designated on logical unit 4 for inclusion, a 
blank data field for KMPR results in the constituent not being included 


in the analysis. A blank card terminates this data type. 


Two cards specifying values for the astronomical arguments SO,HO,PO,ENPO, 


PPO,DS,DH,DP,DNP,DPP in the format (5F13.10). 


SOQ = mean longitude of the moon (cycles) at the reference time origin. 

HO = mean longitude of the sun (cycles) at the reference time origin. 

PO = mean longitude of the lunar perigee (cycles) at the reference time 
origin. 

ENPO = negative of the mean longitude of the ascending node (cycles) at 


the reference time origin. 

PPO = mean longitude of the solar perigee (perihelion) at the reference 
time origin. 

DS,DH,DP,DNP,DPP are their respective rates of change over a 365 day 


period at the reference time origin. 


At least one card for all the main tidal constituents specifying their 
Doodson numbers and phase shift along with as many cards as are necessary 
for the satellite constituents. The first card for each such constituent 
is in the format (6X,A5,1X,613,F5.2,14) and contains the following 


information: 


KON constituent name, 


II,JJ,KK,LL,MM,NN = the six Doodson numbers for KON, 


SEMI the phase correction for KON, 


NJ. the number of satellite constituents. 


A blank card terminates this data type. 


iv) 


If NJ > 0, information on the satellite constituents follows, three 


satellites per card, 


in the format (11X,3(313,F4.2,F7.4,1X,11,1X)). 


For each satellite the values read are: 


LDEL ,MDEL,NDEL = 


PH = 


Be = 


the last three Doodson numbers of the main 

constituent subtracted from the last three Doodson 

numbers of the satellite constituent; 

the phase correction of the satellite constituent 

relative to the phase of the main constituent; 

the amplitude ratio of the satellite tidal potential 

to that of the main constituent; 

1 if the amplitude ratio has to be multiplied by the 
latitude correction factor for diurnal constituents, 

2 if the amplitude ratio has to be multiplied by the 
latitude correction factor for semidiurnal consti- 
tuents, 

otherwise if no correction is required to the amplitude 


ratio. 


One card specifying each of the shallow water constituents and the main 


constituents from which they are derived. The format is 


(6X,A5,11,2X,4(F5.2,A5,5X)) and the respective values read are: 


KON 


NJ = 


COEF , KONCO = 


the name of the shallow water constituent, 
the number of main constituents from which it is 
derived, 


combination number and name of these main constituents. 


The end of these shallow water constituents is denoted by a blank card. 


Logical unit 4 contains six types of data. 


ty One card for the variables IOUT1,RAYOPT,ZOFF,ICHK in the format 
le eras Ceeher lO Os LZ }s 


IOUT] = 6 if the only output desired is a line printer listing 
of results, 
= 2 if both analysis card output and listing are desired. 
RAYOPT = Rayleigh criterion constant value if different from 
sald 
ZOFF = constant to be subtracted from all the hourly heights. 
ICHK = 0 if the hourly height input data is to be checked 


for format errors, 


otherwise if this checking to be waived. 


ii) One card for each possible inference pair. The format is 
(2(4X,A5,E16.10),2F10.3) and the respective values read are: 
KONAN & SIGAN = the name and frequency of the analyzed constituent to 


be used for the inference, 


KONIN & SIGIN = the name and frequency of the inferred constituent, 
R = amplitude ratio of KONIN to KONAN, 
ZETA = the Greenwich phase lag of the inferred constituent 


subtracted from the Greenwich phase lag of the 
analyzed constituent. 


These are terminated by one blank card. 


iii) One card for each shallow water constituent, other than those in the 
standard 69 constituent data package, to be considered for inclusion in 
the analysis. The Rayleigh comparison constituent is also required and 


the additional shallow water constituent must be found in data type i) 


vi) 


of logical unit 8 but have a blank data field where the Rayleigh 
comparison constituent is expected. The format is (6X,A5,4X,A5) and 


a blank card is required at the end. 


One card in the format (11,1X,412,2X,412) specifying the following 
information on the period of the analysis: 
INDY = 8 indicates an analysis is desired for the 


upcoming period; 


0 indicates no further analyses are required; 
1HH1,1DD1,1MM1,1YY1= hour, day, month and year of the beginning of the 
analysis (measured in time ITZONE of input data v)); 


THHL,1DDL,1MML,1YYL= hour, day, month and year of the end of the analysis. 


One card in the format (11,4X,A5,3A6,A4,A3,1X,212,13,12,5X,A5) containing 


the following information on the tidal station: 


INDIC = 1 if J card output is desired (see 1.4), 
= otherwise if not; 
KSTN = tidal station number; 
(NA(J),J=1,4) = tidal station name (22 characters maximum length); 
ITZONE = time zone of the hourly observations; 
LAD,LAM = station latitude in degrees and minutes; 
LOD,LOM = station longitude in degrees and minutes; 
IREF = reference station number. 


The hourly height data cards contain the following information in the 
POV taMbei Ne Oe kao lee One le 
KOLI 


1 or 2 indicates whether this specific card is the 
first or second one for that day, 


= otherwise indicates a non-data card which is ignored. 


JSTN = tidal station number. 

ID, IM,IY = day, month and year of the heights on this card. 

(KARD(J),J=1,12) = the hourly heights in integer form. The final 
constituent amplitudes are in units 1/100 of those 
for the hourly height. Missing values should be 
specified as a blank field or 9999. 


When KOLI=1 the first hourly height on the data card is assumed to be 
at 0100 hr and when KOLI=2 it is assumed to be at 1300 hr. Although it is 
not necessary, for reasons outlined in 2.3.1, it is recommended that these 


observations be recorded in Greenwich mean time. 


After the initial analysis of a computer run is completed, control 
returns to input iv). Successive cards are read then until either a 0 or 8 
value is found for INDY. 


The hourly height data cards need not begin and end so as to include 
exactly the analysis period. The program ignores data outside this range. 
However if more than one analysis is desired from a single job submission and 
hourly height data cards do extend beyond the first analysis period, care should 
be taken to ensure that one of these cards does not have KOLI=0 or blank, 
otherwise the job will be terminated. This is because all successive cards 
after the one containing the last hour of the desired analysis period are read 


in input iv) format. 


1.4 Output 


Four logical units are used for the output of results from the 
tidal heights analysis program. Device number 6 is the line printer; 2 and 7 
are files used for analysis punched card, and alphatext or flexowriter J 
punched card outputs respectively; and 9 is a temporary disc file used for 


storing preliminary punch output and program termination. 6 and 9 are 


required for all program runs whereas the use of 2 and 7 is controlled by 


the input variables IOUT] and INDIC which are read from device 4. 


When IOUT] is 6, INDIC is other than 1, and there are no inferred 
constituents, the only output is two pages on the line printer. The first 
of these lists the constituents included in the least squares fit, their 
frequencies in cycles per hour (although eight decimal places are given, 
depending on computer accuracy, less than this number may be significant), 
the C and S coefficient values (see 2.2.1) measured in units 1/100 of those for the 
hourly heights, and their respective standard deviation estimates. It also 
specifies the number of hourly height observations (excluding gaps) within 
the analysis period, the average and standard deviation of the Original 
observations, the root mean square residual error, and the matrix condition 
number. In the columns titled AL, GL, A, and G, the second page respectively 
lists the amplitudes and phases (degrees) obtained for each constituent from 
the C and S coefficient values, and the same amplitudes and phases after nodal 
modulation and astronomical argument adjustments. The initial and final hour 
of the analysis are also specified along with the Rayleigh criterion constant 
(‘separation'), the midpoint of the analysis period, the total number of 
possible hourly observations in the analysis period, and the total number of 
possible observations used in the analysis. This last value includes gaps in 
the record and is the largest odd number less than or equal to the total 
number of possible hourly observations (if the total number of possible hourly 
observations is an even number, the last hour is ignored). If there is at 
least one inferred constituent, page 2 results are repeated with the inclusion 
of inferred constituents and appropriate adjustments to the constituents from 
which the inferences were made. Appendix 3 lists the final page of results 


obtained from the input values of Appendix 2. 


The only effect of changing the value of IOUT] to 2 (regardless of 
INDIC's value) is to store on file 2, the same information as the second (and 
third) page(s) of the line printer. The list of constituent names, amplitudes 
and Greenwich phase lags begins on line 5 of this file and is in the correct 
format for input to the tidal heights prediction program, namely 


(5X,A5,28X,F8.4,F7.2). 


When INDIC is 1, J card output for the alphatext or flexowriter is 
produced on logical units 6 and 7. This should not be required by the general 
user as it is simply a reorganization of the results in a manner so that they 
can be used for further processing, such as lunitidal computations and the 
production of tide tables. On the line printer these results are given on 
two pages following the standard output. The first of these lists the phase 
and amplitude results in a specific constituent order while the second page 
rearranges them according to increasing frequency and prints them in a format 
Similar to that of device 7. Blank lines separate different constituent 
groups, and on the card punch, these groups are padded if necessary by a 


blank J card to give them an even number of elements. 


1.5 Program Conversion, Modifications, Storage, and Dimension Guidelines 


The source program and constituent data package described in this 


manual were tested on the IBM 370-168 computer installation at the University of 
British Columbia and the UNIVAC 1106 at Institute of Ocean Sciences, Patricia Bay. 


Although as ‘much of the program as possible was written in basic ASCII FORTRAN, 
some changes will have to be made before the program and data package can be 


used on other installations. These may include: 


iii) 


vi) 
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converting the character code from EBCDIC to BCD, 


switching the numeric-character to integer conversion technique in 


subroutine INPUT from DECODE to the local equivalent. 


deleting some or all REAL*8 declaration statements. These are used 

either to permit alphanumeric strings longer than 4 characters (which is 
the ASCII FORTRAN single precision word length) to be read (e.g. constitu- 
ent names are read in A5 format), or to gain computational precision in 

the average, standard deviation, and sine/cosine iteration formula 
calculations of subroutine SCFIT2. CDC installations should not require 
double precision variables for either of these cases because their single 
precision is about 1.5 that of IBM and UNIVAC, and because there are no 
alphanumeric strings longer than 6 characters (which is the CDC word 


length) in the program. 


altering the variable list structure for ENTRY statements and reference 


to them. 


changing some or all of the logical unit (file reference) numbers from 
their present values (namely 2,4,6,7,8,9) in order to conform with 


local machine restrictions. 


deleting, or changing the form of, the END FILE statement in subroutine 


DGM388 which writes an end of file mark before rewinding device 9. 
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The program in its present form requires approximately 5000, and 
52000 UNIVAC words, for the storage of its instructions, and arrays 
respectively. A large part of this is due to AS, the 170 by 170 matrix 
resulting from the least squares fit for constituent amplitudes and phases 
(see 2.2.1). If memory requirements are restrictive on a particular 
installation, array storage can be cut to approximately 38000 words by 
storing only the upper triangle of this matrix in a one dimensional array. 
A program version that does this is available upon request, however it will 
have somewhat slower processing times than the one described in this manual 


(see 2.2.2 for details). 


Changing the number or type of constituents in the standard data 
package may require some alterations to the analysis program. The arrays 
NAME and ARRAY in subroutines SETTAB and DGM388 respectively, are used for 
producing J card output, and their contents to have been designed 
specifically for the present data package. If J card output is desired, 
appropriate changes to these arrays should accompany any to the data 
package. If constituents are added to the standard data package, the 
dimensions of several arrays may have to be altered. Restrictions on the 


minimum dimension of such arrays are now given. 


Te 


Let 

MTOT be the total number of possible constituents contained 
in the data package (presently 146), 

M be the number of constituents considered for inclusion 
in the analysis (presently 69 + the number of shallow 
water constituents specifically designated for inclusion), 

MCON be the number of main constituents in the standard data 
package (presently 45), 

MSAT be the sum of the total number of satellites for these 
main constituents and the number of main constituents 
with no satellites (presently 162+8), 

and MSHAL be the sum for all shallow water constituents, of the 


number of main constituents from which each is derived 


(presently 251). 


Then in the main program, arrays KONTAB, FREQ, and KMPR should have 
minimum dimension MTOT; arrays KON,C,S,SIG,ERC,ERS,A,EPS,KO,AA and GD should 
have minimum dimension M; array AS should have minimum dimension 2M-1 by 2M-1; 
array NKON should have dimension at least as large as the number of extra 
shallow water constituents specifically designated for analysis inclusion (its 
present maximum is 15); and array Z should be large enough to contain the 


hourly heights (and gaps) in the analysis period (its present maximum is 


375 days). 


In subroutine INPUT, because arrays Z,AS,A, and EPS are in a common 


block, they should be dimensioned the same as in the main program. 


ihe 


In subroutine VUF, arrays KON,VU,F and NJ should have minimum 
dimension MTOT; arrays II,JJ,KK,LL,MM,NN, and SEMI should have minimum dimension 
MCON; arrays EE,LDEL,MDEL, NDEL, IR, and PH should have minimum dimension MSAT; 
and KONCO, COEFF should have minimum dimension MSHAL. 


In subroutine INFER, arrays KONAN ,KONIN,SIGAN,SIGIN,R, and ZETA can 


presently accommodate a maximum of 10 inferred constituents. 


In subroutine SETTAB, arrays KO,AA,GD and NKON should be dimensioned 


as in the main program, while NAME should have minimum dimension M. 


In subroutine DGM388,ARRAY should have minimum dimension M; and 
arrays IZONE,CT,DATAN,NUM,JJ,STATNO,G, and H should have minimum dimensions 
M+1. 


In subroutine SCFIT2, arrays X,A,AM, and EPM are in COMMON with 
Z,AS,A, and EPS of the main program and hence should have the same dimensions. 
Arrays CW,SW,RHSC, and RHSS should have minimum dimension M, and array RHS 
should have minimum dimension 2M-1. AC and AS should have the same first 
dimension as A and care should be taken that through their equivalence rela- 


tionships, neither AC and AS, nor RHSC and RHSS overlap. 


Finally, in subroutine CHLSKY, arrays A and F should have minimum 


dimensions 2M-1 by 2M-1,and 2M-1 respectively. 
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2: TIDAL HEIGHTS ANALYSIS PROGRAM DETAILS 


2.1 The Constituent Data Package 


2.1.1 The astronomical variables 


The astronomical variables required by the tidal analysis program were 
used by Doodson [1921] in his development of the tidal potential. From them one 
can calculate the position of the sun or moon, and hence the tide generating 


forces, at any time. These variables are: 


S(t) = mean longitude of the moon, 

H(t) = mean longitude of the sun, 

P(t) = mean longitude of the lunar perigee, 

N'(t) = negative of the longitude of the mean ascending node, 
and P'(t) = mean longitude of the solar perigee (perihelion). 


For H, N' and P' these longitudes are measured along the ecliptic eastward from 
the mean vernal equinox position at time t; while for S and P they are measured 
in the ecliptic eastward from the mean vernal equinox position at time t to the 
mean ascending mode of the lunar orbit, and then along this orbit. Together 

with the rates of change of these variables, t the local mean lunar time, and 

the Doodson numbers for each tidal constituent, one can calculate the constituent 
frequencies, their astronomical argument phase angles V, and their nodal modula- 


tion phase, u, and amplitude, f, corrections. 


The values of the astronomical variables in the constituent data 
package were calculated from power series expansion formulae in the Explanatory 
Supplement to the Astronomical Ephemeris and the American Ephemeris and Nautical 
Almanac [1961], and the Astronomical Ephemeris [1974]. These formulae were 
derived from Newcomb's Tables of the Sun and a revision of Brown's lunar theory 


(used in the development of his Tables of Motion of the Moon) so that it is in 
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accord with Newcomb's. Calculations were done in double precision on an IBM 


computer and have approximately 14 significant digits accuracy. 


In the program, the astronomical variables used in the calculation of 
V(t), u(t) and f(t), are evaluated to the first term in their Taylor expansion; 
e.g. S(t) = S(t.) + (t-t,) Bt) Provided they are consistent with other 
time dependent Ayer eat aatigg and ty can be chosen arbitrarily. However, in 
order to facilitate computations (see 2.2.1), t is chosen to be the central 
hour of the desired analysis, and in order to gain precision to: the reference 
time origin, is taken to be 000 ET! January 1, 1976. This latter date, it was 
felt, would be closer to the analysis period of most records than the previous 


choice of 000 ET January 1, 1901, and hence would yield more accurate results 


via the linear approximation. 


In keeping with the choice of reference time origin and astronomical 
variable specifications, t should be measured in Ephemeris time. However, the 
correction from Universal time is irregular and in most cases small, so it 
has been assumed for computational purposes that all observations are recorded 


ci eet eae es 


The frequencies specified for all constituents have been calculated 
from the rates of change of the astronomical variables and the respective 
Doodson numbers. This means that they too are calculated at the time origin 
January 1, 1976; although there would be very little change were they taken 
at January 1, 1901 instead. | 


‘ephemeris Time (ET) is the uniform measure of time defined by the laws of 


dynamics and determined in principle from the orbital motion of the Earth 
as represented by Newcomb's Tables of the Sun. Universal or Greenwich Mean 
Time is defined by the rotational motion of the Earth and is not rigorously 
uniform. 
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2.1.2 Choice of constituents and Rayleigh comparison pairs 


There is a maximum of 146 possible tidal constituents that can be 
included in the tidal analysis, 45 of these are astronomical in origin 
(main constituents) while the remaining 101 are shallow water conceituents. i 
Because computation time (and cost) of the computer program increases 
approximately as the square of the number of constituents included in the 
analysis, and because for many tidal stations, most of the shallow water 
constituents are insignificant, a smaller standard package was seen as 
adequate for general use. Based on the suggestions of G. Godin, it was 
decided that this package contain all the main constituents and 24 of the 


shallow water. However, provision was made so that other shallow water 


constituents among the 77 remaining could be included if desired. 


The Rayleigh comparison constituent is used for the purpose of 
deciding whether or not a specific constituent should be included in the 


analysis. If Fy is the frequency of such a constituent, Fy is the frequency 
of its Rayleigh comparison constituent and T is the time span of the 


proposed record to be analyzed, then the constituent will be included in 
the analysis only if |Fo - F,| T > RAY. RAY is commonly given the value | 


although it can be specified differently in the program. 


] : A : : " ; 
( iy criterion for selecting these main constituents was to include 
all the diurnal and semidiurnal constituents with Cartwright [1973] 
tidal potential amplitudes greater than 0.00250, along with M3 and 
the most important low frequency constituents. Section 2.1.3 gives 
the analogous shallow water constituent criterion. 


iv, 


In order to determine the set of Rayleigh comparison pairs, it is 


important to consider, within a given constituent group (e.g. diurnal or 


semi-diurnal), the order of constituent inclusion in the analysis as T 


(the time span of the record to be analyzed) increases. Assuming this 


point of view, the specific objectives used when constructing the set 


listed in Appendix 1] were: 


ji) 


iii) 


within each constituent group, when possible, have the order of 
constituent selection correspond with decreasing magnitude of tidal 


potential amplitude (as calculated by Cartwright [1973]), 


when possible, compare a candidate constituent with whichever of 
the neighbouring, already selected constituents, that is nearest 


in frequency, 


when there are two neighbourina constituents of relatively equal 
tidal potential amplitude, rather than waiting until the record 
length is sufficient to permit the selection of both at the same 
time (i.e. by comparing them to each other), choose a representa- 
tive of the pair whose inclusion will be as early as possible. 
This will give information sooner about that frequency range, and 
via inference, still enable some information to be obtained on 


both constituents. 
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The Rayleigh comparison pairs chosen for the low frequency, diurnal, semi- 


diurnal and terdiurnal constituent groups are given in Tables 1, 2, 3 and 
4 respectively. Figures given for the length of record required for 
constituent inclusion assume a Rayleigh criterion constant value (input 


variable RAYOPT) of 1.0. 


201 and SIG] provide an example of objective iii). Because 2Q1 
has a greater frequency separation from Q1 and hence would appear in an 
analysis of shorter record length than SIG1, it was chosen as the 


representative. 


However, it can be seen in several cases, that it was not possible 
or feasible to adhere to all the objectives just outlined. Choosing a 
Rayleigh comparison constituent from the list of those constituents already 
included in the analysis proved to be difficult near the frequency edges 
of constituent groups. Upward arrows indicate failure to uphold this 
objective. 001 is such a case. For it, the potential comparison pairs 
were SO1, Kl and Jl. The first of these would result in both SO1 and 001 
appearing at the same later time than had Jl or Kl been chosen. Hence, 
information about 001 would be unnecessarily delayed. Although, due to 
the tidal potential amplitude of Jl, objective i) is violated with both 
the second and third choices, it was felt that the third was a better 


compromise. With it, 001 only appears 11 hours sooner than Jl. 
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Table 
Shallow Water Constituents in 


Length of 
Record (hr) 
Required 
for 
Constituent 


Shallow 
Water 


Component 
Lengths 
1 for Their 


(hr) of 


5) 


the Standard Data Package 


Main Constituents 
Record 


the 


and 
Required 


Inclusion in Analysis 


24 


K2 is an example of an unavoidable violation of objective i). Because 
jt is so close in frequency to S2, its importance as a major semidiurnal consti- 


tuent does not insure it an early inclusion in the analysis package. 


Because shallow water constituents do not have a tidal potential ampli- 
tude, objective i) does not apply to them. However, based on his experience, 
Godin was able to suggest a hierarchy of their relative importance. A further 
criteria used when selecting comparison pairs for them was that no shallow water 
constituent should appear in an analysis before all the main constituents, from 
which it is derived, have also been selected. Table 5 shows that this has been 
upheld for all the shallow water constituents in the standard 69 constituent 


data package. 


I suggest that the objectives outlined here be employed when choosing 
the Rayleigh comparison constituent for any additions to the list of possible 


constituents to be included in the analysis. 


2.1.3 Satellite constituents and nodal modulation 


Doodson's [1921] development of the tidal potential contains a very 
large number of constituents. Due to the great length of record required for 
their separation, several of these can be considered, for all intents and 
purposes, unanalysable. The standard approach to this problem is to form 
clusters consisting of all constituents with the same first three Doodson 
numbers. The major contributor in terms of tidal potential amplitude lends its 


name to the cluster and the lesser constituents are called satellites. 


The method of analysis uses this main and satellite constituent approach 
in the following manner. The Rayleigh criteria is applied to the main constituent 


frequencies to determine whether or not they are to be included in the analysis. 


<a 


For each of those so chosen, we analyse at its frequency and obtain an apparent 
amplitude and phase. However, because these results are actually due to the 
cumulative effect of all the constituents in that cluster, an adjustment is 

made so that only the contribution due to the main constituent is found. This 


adjustment is called the nodal modulation. 


In order to make the nodal modulation correction to the amplitude and 
phase of a main constituent, it is necessary to know the relative amplitudes 
and phases of the satellites. As is commonly done, it is assumed in this program 
that the same relationship as is found with the equilibrium tide (tidal potential), 
holds with the actual tide. That is, the tidal potential amplitude ratio of a 
satellite to its main constituent is assumed to be equal to the corresponding 


tidal heights amplitude ratio, and the difference in tidal potential phase equals 


the difference in tidal height phase. 


The source of the tidal potential amplitude ratio, as found in the 
constituent data package of Appendix 1, is Cartwright [1971,1973}. Using new 
computation methods and the latest values for the astronomical constants, he 
obtains more accurate results than those from the previously used Doodson 
computations. It should be noted that in several cases (whenever the satellite 
arises via a third order term), this version of the constituent data package 
requires that the amplitude ratio be multiplied by a latitude correction factor. 
This represents a refinement of the previous version where a latitude of 50°N 
was assumed and the multiplication factor was incorporated directly into the 


ratio. 


Phase differences between satellites and main constituents arise when 
the tidal potential development yields different trigonometric terms for these 


constituents. The common convention is to express all terms in cosine form and 
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so an extra -% cycle phase shift is introduced if the term was originally a 
sine. Satellites requiring such a shift are called third order. A further 
% cycle change is also introduced when all negative amplitudes are made 


positive. 


Because several test analyses indicate less consistent results when 
third order satellites are included in the N2 and L2 nodal modulation, Godin 
has decided to delete these from the present standard constituent data package. 
Instead he suggests that the results of analyses with this package should be 
compared with those of previous analyses in order to find the most suitable 


adjustment for these constituents. 


The only other main constituents that do not have all their satellites 
included for nodal modulation are the slow frequency constituents. For them, 
no satellites are specified. Because low frequency noise may be as much as an 
order of magnitude greater than the satellite contributions, and Mm, MSf and 
Mf when they are detectable are often of shallow water origin, the effect of 
making corrections for the expected satellites would be to obscure further, 


rather than clarify the actual low frequency periodic signal. 


Section 2.3.2 gives further details on the nodal modulation correction. 


2.1.4 Shallow water constituents 


Shallow water tidal constituents arise from the distortion of main 
constituent tidal oscillations in shallow water. Because the speed of propaga- 
tion of a progressive wave is approximately proportional to the square root of 
the depth of water in which it is travelling, shallow water has the effect of 
retarding the trough of a wave more than the crest. This distorts the original 
Sinusoidal wave shape and introduces harmonic signals that are not predicted in 


tidal potential development. The frequencies of these derived harmonics can 


27 


be found by calculating the effect of non-linear terms in the hydrodynamic 
equations of motion on a signal due to one or more main constituents (see Godin, 


[1972] pages 154-164 for further details). 


The shallow water constituents chosen for inclusion in the standard 
69 constituent data package were suggested by G. Godin. They are listed in 
Table 5 and are derived only from the largest main constituents, namely 
M2,S2,N2,K2,K1 and 01, using the lowest types of possible interaction. The 
77 additional shallow water constituents that can be included in the analysis 
if so desired are derived from lesser main constituents and higher types of 
interaction. In the constituent data package listing of Appendix 1, they can 


be spotted by their lack of a Rayleigh comparison constituent. 


When shallow water effects are noticeable, main constituents, if they 
are close in frequency, may coexist or be masked by constituents of non-linear 
origin. The resultant nodal modulation will be due to the pair and thus will 
not coincide to the calculated modulation of the main constituent. In 
suspected cases, the effectiveness of nodal corrections in a series of 
successive analyses will indicate the presence of pairs or emphasize the pre- 
dominance of one constituent over the other. The following table (taken from 
unpublished notes of Godin) lists compound constituents which may coexist 
with or mask constituents of direct astronomical origin. In all cases except 
SO1 and M03, the main rather than the compound constituent is included in the 


standard constituent data package. 
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TABLE 6 


Shallow Water Constituents that may Mask Main Constituents 


Compound Constituent 
which may coexist at or 
near its frequency 


Main Constituent 


Q] 


01 MK1 * 
TAU MP1 * 
Din Semen ws io + 
PI SKI * 
K1 MO1 
J MQq 


* The modulation or frequency of the compound constituent is 
sufficiently different that the pair could be separated if 
a long enough record of high precision were available. 


29 
2.2 The Least Squares Method of Analysis 


2.2.1 Formulation of the problem 


The first stage in the actual analysis of tidal records is the least 
Squares fit for constituent amplitude and phase. If the tidal record is 
of minimum length 13 hours, the present program and data package insure that 
the constant constituent ZO is always included in the analysis. If 
dj, £0rj a b3M are .the frequencies (cycles/hr) of the other tidal constituents 
chosen for inclusion in the analysis by the Rayleigh criterion, then the 
problem is to find the amplitudes A. and phases dj of the function 
Ce.* dy A. cos (2n(ost, - o3)) that best fit the series of observations 
Assuming N > 2M + 1 we see that it is impossible to solve 
the system y(t; )=C, - eS A. cos (2n(o.t, - 03) exactly because it is over- 
determined. Hence, fee necessary to adopt a criterion which will enable 
unique optimum values for the parameters A. and os to be found. The most 
common optimization criterion used, and the one chosen here, is the least 
Squares technique. 


Re-expressing a 


fn A. cos (2n(o.t. - 6;)) as 


ie (C, cos (210 .t,) + S; sin (200 5t;)), 


where As = vc? + S57 and = arctan Sy lls so that the fitting function 
is linear in the parameters oy and C. and hence more easily solved; and 
rewriting y(t.) as Ya3 the objective of the least squares technique is to 


minimize 


(1) 


In order to minimize the loss of accuracy due to round off, the average 
of the hourly heights observations is subtracted from all original values. 
The y(t.) values mentioned in all computations henceforth are actually 
the resultant deviations. At the end of all calculations, Cy is adjusted 
by this mean value. 
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M in 2 4 
-Cy- 2 (C; cos 2ma,t; +S; sin 2no;t;)) 


a SS. j=1,M. 
for Cy and all C595; a 


This is done by solving the following 2M + 1 simultaneous equations for 


J=1,M 

px io EY ve (y, =, " C, cos 2no.t, - a S, sin 2no,t;) (-1) 3 

0-2 = 2") (y, - C5 - NE Cj cos anayt, - MY S, sin 2noyt,) (- cos 2na,t,) 5 
5 i=] j=l j=l 

Gs= 38; Sy ” (Yrer'C Sie te C, Cos ano st J ie S; sin 2no sts) (- sin ano jt). 


This results in the matrix equation Bx = y of Figure 1 (page 33). 


Gaps in the data record (i.e. missing hourly observations) are easily 
handled by the least squares method because it is not necessary that the 
observation times t. i = 1,N be evenly spaced. For example, if the analysis 
covers the total time period of 100 hr but hours 50 to 74 inclusive are 
missing, then tsg will correspond to the 75th hour. However, because the 
following identities which simplify the summations require that the observ- 
ation times be evenly spaced, it is necessary that each of the matrix terms 
be calculated as the sum of contributions over the data periods that contain 
no gaps. Assuming that [ng,n,] is the hour range of a section of record 
containing no gaps, we can substitute ty = k in the matrix coefficients 


expressions since the times are at successive hours. 


Using the relationships 
cos a cos b = 4(cos(atb) + cos(a-b)) 


Sin a sin b = &(cos(a-b) - cos(atb)) 


Sin a cos b = }(sin(atb) + sin(a-b)), 


Si 


the formula for the sum of a geometric series, namely 


n n+] 
a+ ar t expo ane s.alr 


r" 1)/(r if 1 


3 é d : 1x 
and expressing cos x and sin x as the real and imaginary parts of e 


, we 


obtain the identities: 


OL). EOS aK Xaeke IT) ea ol cos [(™*%0)x| sin Capa 
k=Nno 2 2 


and 1) sin kx = sin (aa sin [(M*"0)x|/sin (x/2) 
k=NnNo 2 2 


Hence, the summation expressions in the least squares matrix can be simplified 


(with regard to computer execution time) as follows. 


may cos (2n0,k) cos (2no2k) = 1/2 ny [cos (2mk(o,+0,)) + cos (2rk(s}-e2))| 
k=Nng k=no 


1/2 (sin ((ny-ngt])m (oy+02)) cos ((nytng)t (o,+0,)) / sin t (a, 40) 


+ 


sin ((ny-ngt1)m (cy-05)) cos ((nytng)t (0y-02)) / sin © (a1-32)) 


My Sin (2n0,k) sin (2n05k) = 1/2 my [cos (2nk(o,-05)) - cas (2nk(s2*22))| 
k=Ng K=Nno 


= 1/2 (sin ((ny-not1)m (a;-a2)) Cos ((nytng)m (oy-2)) / sin w (33-09) 


BiSin ((ny-ngt1) 1 (oy to5)) sin ((ny +N) (o,t0>)) LST e: (or ee oye 


ny sees 
,: Sin(2iro,k) cos (21o,k) = ehhZ any [sin (2rk(c,+0,)) - Sn ek s -3,))| 
k=no K=Nng | : 


= 1/2 (sin ((n)-n +1 )« (o,40,)) sin (nj +n, )n (gta) etesin « (3, +9,) 


+ sin ((n,-n,+1) 1 Lorre juss ({n, +n.) (o,-0,)) / sin + (c,-<,)) 
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With these substitutions made in Fig. 1 we have the least squares 


matrix equation Bx = y generated in subrouting SCFIT2. 


Because B is symmetric it is sufficient to store only its upper 
triangle consisting of 2M* + 3M + 1 elements instead of the entire matrix of 
(2M + 1)2 elements. However this saving does not come without its price as 
the complicated subscripting scheme required to keep track of these elements 
when stored in a one dimensional array significantly affects the solution 
computation time (see 2.2.2 for details). Assuming that this program will 
not have storage problems on most computer installations, it was decided to 


waive the storage saving in favour of improved computational time. 


Partitioning the matrix equation Bx = y into the form 
ie 
Ysto)> 
where Bi,, By2, Boy, Boos Gs Ss Yes Ys have dimensions (M+1) KX (M+1), 
(M1) XMS MK OUMET) MX MS (MET Ket Sm XOr MeL) Xe Too Mex TV respectively. 


it is easily seen when no=-n,, that B,> and Bo, become zero matrices and 


Bia Bio 


lo 
il 


Ys) 


Bo, Boo 


two smaller matrix equations, B,,c = y, and Boos = y, result. The combined 
computation time to solve these equations is less than that of the original 

(see 2.2.2) so it is desirable to attain this condition when possible. Since 
the time origin of the hourly observations is arbitrary provided it is 
consistent with that of the astronomical argument V, we can attain the 

desired condition for a record with no gaps by choosing the central hour of 

the record as the origin. (This requires that the total number of observations 
be odd and is satisfied by ignoring the last observation, if the total is 

even). Although there is generally no corresponding matrix simplification 

in the case of a record with gaps, for consistency with the foregoing choice, it 


is convenient to choose the central hour of the record universally as the time 


origin. 


of 


FIGURE 1 
The Matrix Equation Bx = y Resulting from the 


Least Squares Fit for Constituent Amplitudes and Phases 
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2.2.2 Solution of the matrix equation, the condition number, 
and statistical properties 


Most of the discussion and development of the Cholesky factorization 
algorithm introduced in this section is taken directly from Forsythe and Moler 
[1967]. Although all results and discussion are now stated only for the 
matrix B and the equation Bx = y, they apply as well for the partitioned systems 


BYd> B,C = Ye and Boos BooS = Yo. 


In addition to symmetry, a useful property of the matrix B in Fig. 1, 
is that it qualifies the matrix equation for a simpler solution than were it a 
general matrix, is its positive definiteness. This property requires that 


for all (2M+1) X 1 dimensional vectors x # 0, x!Bx > Ge 


The positive definiteness of B can be demonstrated by considering the 
overdetermined matrix equation y = Ax + e resulting from the system of 
equations y(t.) op Gm i (C. 
the vector x! = (GPO). PHB StiSs seen Seis y Men (yl tr) a ene aya and 


cos eno st. + S; Sin eno .t.) rie for i = 1,N where 


e is the vector of residuals. It is easily seen that Ala = B, and so for any 


x # 0, 


.2, where xlal = z! =f by a ee 


x'Bx = x AlAx = z!7 = by z 
ae is Mae ibe 


Nn) 

It is worth mentioning that the overdetermined system y = Ax +e can 
be solved in many ways, depending on the criterion chosen for minimizing e. 
For our purposes, those methods which solve the system without changing 
the form of the matrix are impractical from a Storage, processing time and 
rounding error point of view because the first dimension of A(= the number of 


hourly observations) is commonly 9000. However, minimizing ele is equivalent 


to the least squares criterion adopted here. 
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An important result for any positive definite symmetric matrix B is 
that it can be uniquely decomposed in the form B = ene where G is a lower 
triangular matrix with positive diagonal elements (1). Expanding this 


relationship leads to the matrix element equalities: 


oe Jy marie 
id alencnecalk 


sid eo 
bs, ey IiKI5k for alvey Sy? 


The algorithm resulting from using these equations in the proper order to find 
the elements of G is known as Cholesky's square root method for factoring a 
positive definite matrix (also attributed to Banachiewicz; see Faddeev and 
Faddeeva [1963]). Unlike other matrix decomposition methods such as Gaussian 
elimination, it does not have to search for and divide by pivots. Such 
techniques must insure that the reduced matrix elements are not too large so 
that rounding errors and loss of accuracy do not occur. In Cholesky's method 
however, we can see that 1955 < Yb;, for all i,j and so upper bounds for the 


elements of G always exist. 


Prior to revision of SCFIT2 the constituent ZO was treated like all 
others in the sense that both Cg and Sy were found. In order to insure that 
ZO's phase was zero, Sg was forced to zero by setting the corresponding 
diagonal element of B to a very large number. Although this did not affect 
the symmetry or positive definiteness of the matrix, there was a loss of 
accuracy (which is more evident on computers, such as IBM, with a smaller 


number of single precision significant digits) due to the corresponding 


(1) 


If B is symmetric but not positive definite a similar decomposition exists. 
However some elements of G may be complex or, in the degenerate case, zero 
along the diagonal. 
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large element in the reduced matrix. For this reason and the slightly 
reduced array storage required, in the revised subroutine, only Co is now 


sought in the least squares fit. 


Once B has been decomposed into the upper and lower triangular 
matrices, it is a relatively easy matter to solve the matrix solution. This 
is done by breaking down the equation 6g!x = y into Gb = y and gly anbe 
Because of the triangular nature of G, these equations can be solved by 


forward and backward substitution for b and x respectively. 


The amount of arithmetic in a matrix algorithm is usually measured 
by the number of multiplicative operations (i.e. multiplications and divisions) 
used, since there are normally approximately the same number of additive 
operations. For a matrix of dimension n x n, the Cholesky factorization 
algorithm requires n square roots and approximately 1/6 n° multiplications. 
This compares favourably with the 1/3 n? multiplications required by Gaussian 


elimination (Wilkinson [1967]) to produce a triangular matrix. 


Wilkinson [1967] suggests a factorization of B into LoL', where Lams 
a lower triangular matrix and D is a positive diagonal matrix, that involves 
no more multiplications than Cholesky and avoids the square roots. However, 
assuming that the time ratio of a square root operation to a multiplication 
is 15:1 (approximate ratio for IBM 370-168) and that all 69 constituents in 
the data package are included in the analysis (i.e. n = 137) the time saved 
by eliminating the square roots is only 0.5%. Furthermore some of this gain 
would be replaced by time required for storing and retrieving information from 
the additional matrix D, and for the n additional division operations each 
time a solution is calculated by forward and backward substitution. Hence, 


the factorization was not adopted in the present program. 


oy 


Because the time required for the factorization of B varies as the 
cube of the number of unknowns, an approximate four-fold time reduction 
Should result when the tidal record has no gaps and the partitioned rather 
than the original matrix equations are solved. However as the following 
table of execution times for sections of subroutine SCFIT2 demonstrates, 
Significant improvements can also be cited in the time required for 
matrix generation, and error calculation. The values shown were obtained 
on an IBM 370-168 computer with a 34 constituent analysis of a 38 day tidal 


record. 
TABLE 7 
Comparison of Processing Times between the Partitioned 
and Non-partitioned Matrix Equation Solutions 


Partitioned matrix : Non-partitioned 
system times (sec.) | matrix system (sec. ) 


itunes #2 ted! gonad Gy | 
Parameter initializations | 
& right hand generation 347 | . 346 
Matrix generation 059 | 2178 
Matrix factorization 049 | . 146 
Solution | .010 | .018 
Error calculations | 5 Was: -403 


The advantage of storing the lower triangle of B in a two dimensional 
array and thus avoiding the complicated subscripting scheme required for one 
dimensional storage is illustrated in TABLE 3. The values shown were obtained 
from solving a matrix equation of order 110 resulting from a yearly tidal 
analysis. For comparison, the times required when using the Forsythe and 
Moler technique (Gaussian elimination with partial pivoting, forward and back- 


ward substitution) with two dimensional matrix storage are also included. 
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TABLE 8 


Comparison of Processing Times for the Cholesky and Forsythe & Moler Methods 
of Solving the Least Squares Matrix Equation 


| PROCESSING TIMES (SEC.) 


ele method | Cholesky method} Forsythe and 


| with two with one Moler method 
' dimensional dimensional with two 
' array storage array storage dimensional 


| 
| 
| array storage 


Se eet ett OE aie ee ae oat OG ee Ss DR es Bete ee 
Matrix factorization | wes | wena Laeg 2.38 
Solution ! 053 059 .047 


A rough indication of the roundoff difficulties associated with 
solving the equation Bx = y is given by the matrix condition number. Although 
several different definitions for a condition number exist, an appropriate 
one for our purposes, in the sense that it pertains to least squares matrices 
and is easily calculated, is specified by Davis and Rabinowitz [1961]. Its 


development is as follows. 


ECV AID pee, ote b,} are n-dimensional vectors such that the matrix 
by Chameziee Dat | bb ras bib, 
ae 1 ' 
pee = | 
I 1 i 
' ! ' 
b. bibl ...-. b,-b, ; 


then it can be shown that 0 < det(B) < lb, | Ib, || 1* 4 ical where if 

b. = (bs. a. bin) the norm i = sana Furthermore det(B) = 0 
i) 

if and only if the vectors are linearly dependent, and det(B) = Ib. | ad TY 


Hb, || if and only if they are orthogonal (i.e. b,-b. =e Olfor Wa P.us This 


determinant is known as the Gram determinant of the SVS tem 1/056 tans bi} and 


is the square of the n-dimensional volume of the parallelepiped whose edges 


are these vectors. 
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Since it can be shown that all least squares matrices can be expressed 
in this manner, this result can be applied to our situation. In particular 
when the vectors are normalized so that Ib.1] = 1, the actual value of det(B) 
will always be bounded and provide a measure of the linear independence of the 
system, and hence roundoff difficulties encountered in solving the equation. 

A value close to 1 will mean near orthogonality, a virtually diagonal matrix 
for B, and thus an easy solution. On the other hand, a value close to 0 will 
mean that at least two rows are near scalar multiples of one another, and thus 
greater accuracy problems will occur when their difference is calculated 
during the equation solution. 


For our particular case observe that det(B) = det (GG! } = (detG)? = 
n 


II . 
2 _Sii, and that B can be written as 


Sy Fc Gi Gy 
gg! = 
Jn Sn Jn Zn 
911 mee Snl | 
where g! x 922 Sno Efe re 92 9) 
| Or aie caste 0g, 
Since b.. = Jy gh, Hg .{| = vb.. and the determinant of the matrix resulting 
AP k=] jk 3] dl 


from normalizing the ae Vectors 15 r ne a ibiaae The square root of this value 


“i 


is the volume of the n-dimensional pa mre mia 1s whose edges are these 
normalized vectors and is the quantity calculated as the condition number of 


the matrix B. 
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The statistical properties of the least squares fit solution can be 
found in any analysis of variance or regression model text. They are outlined 


briefly as follows. 


Reverting to the overdetermined problem statement, the least squares 


t 


objective can be stated as finding the vector x in y = Ax + @ such that ee 


X 
is minimized. This yields the solution x = (A‘A)~} aly, 

The total sum of squares is y!y and the sum of squares due to 
regression is xlaly. Their difference is the residual error sum of squares 
and this difference divided by the degrees of freedom in the fit is the 
residual mean square error (MSE). ‘Degrees of freedom' is the difference 
between the number of hourly observations (excluding gaps) and A the number 
of parameters fit in the analysis. If there were M constituents including ZO 


chosen for the analysis,the degrees of freedom would be N - 2M + 1. 


If it is assumed, as is commonly done, that the vector e is distri- 
buted normally with 0 standard deviation and o“I variance, where I is the 
unit diagonal matrix, then the variance of x is (Aare o*. Since the mean 
Square residual error is an unbiased estimator for o%, an estimate of the 


standard deviation of x,, the ith element of x, is 


1/2 


: Tyy-1 

((u,! (ATA)"? u.) MSE) 
where uy; is the vector with 1 in the ith position and 0's elsewhere. 
2.3 Modifications to the Least Squares Analysis Results 


2.3.1 Astronomical argument and Greenwich Phase Lag 


Instead of regarding each tidal constituent as the result of some 


particular component of the tidal potential, an artificial causal agent can be 


4) 


attributed to each constituent in the form of a fictitious star which travels 
around the equator with an angular speed equal to that of its corresponding 
constituent. Making use of this conceptual aid, the astronomical argument, 
V(L,t), of a tidal constituent can then be viewed as the angular position of 
this fictitious star relative to longitude L, and at time t. Although the 
longitudinal dependence is easily calculated, for historical reasons L is 
generally assumed to be the Greenwich meridian, and V is reduced to a function 


of one variable. 


The Greenwich phase lag, g, is the difference between this astrono- 
mical argument for Greenwich and the phase of the observed constituent signal. 
Its value is dependent upon the time zone in which the hourly heights of the 
record were taken. This means that when phases at various stations, not 
necessarily in the same time zone, are compared, they must be reduced to a 
common zone in order to avoid spurious differences due to difference relative 
times. Specifically, if o is the constituent frequency and g(j+a,) and g(j) 
are the Greenwich phase lags evaluated for time zones Jt, and j respectively 


(e.g. Pacific Standard Time is +8), then 


g(j+a5) = 9(5) + (d,)o. 
Although these adjustments are easily calculated, they can be tedious because 
each constituent must be handled individually. Therefore, to avoid possible 
misinterpretation of phases from nearby stations or subsequent phase 


alterations, it is suggested that all observations be recorded in, or 
converted to, GMT. 


The calculation of g (see 2.3.3) requires that the astronomical 
argument need only be evaluated at one time, the central hour of the analysis 


period. For a particular main constituent, it is calculated as 


Va dig pt oes ako Soha Rech man ot ngP 
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where igsdjosko.%osMg.Nqo are the Doodson numbers of the constituent and 
S,H,P,N',P' are the astronomical variables defined in 2.1.1. 1, the number 
of mean lunar days from the absolute time origin of 000 ET January 1, 1976, 
is calculated as sum of the local mean solar time from this origin and 


(H-S), and so need not be read from data cards. 


For shallow water constituents, the astronomical argument is calcu- 
lated as the linear combination of the coefficient number and the astrono- 
mical argument of the main constituents from which it is derived. For example, 


Ve,e7w - V,. and V 


Vusn2 = “m2 + “so > "no OMK5 = 2 Yo + Vyq- 


2.3.2 Nodal corrections 


Most of this section has been taken from the unpublished notes of 
G. Godin which were written subsequent to Cartwright's [1971,1973] recalcu- 
lation of the tide generating potential. The material presented here is 


intended to give greater detail than that of section 2.1.3. 


Due to the presence of satellites to the major contributor in a 
given cluster, it is known from tidal potential theory that the analyzed 
Signal found at the frequency G3 of the main constituent is actually the 


result of 


a; sin (V - 9) + } Axe az, Sin (Viz - Diy) + ) Ae Az, COS (Vio - Ii) 


for the diurnal and terdiurnal constituents of direct gravitational origin, and 


- g.,) 


a; cos (V, = 9;) +.) Aik as, COS (Vey - Jy) + ) a ai, sin (V “i 


k 


for the slow and semidiurnal constituents. a, g and V are the true amplitude, 


je 


Greenwich phase, and astronomical argument at the central time of the record 
for all the constituents. Single j subscripts refer to the major contributor 


while jk and jx subscripts refer to satellites originating from tidal potential 
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terms of the second and third order respectively. A is the element of the 


interaction matrix resulting from the interference of a satellite with the 


main constituent. 


It is the convention in tides and an assumption for our least squares 


Fit that all constituents arise through a cosine term and positive amplitude. 


That is, the contribution for a constituent whose astronomical argument is Le 


and whose Greenwich phase is 95> is expected to be in the form a COS) (Vi 


7OF 2. 
or a; 


= 


57 95) 
0. However, the diurnal and terdiurnal constituents, assuming that 


they are due to second order terms in the tidal potential, actually arise 


through a b, sin (V. 


Pas 9.) term where b. may be negative. Hence, a phase 


correction (variable SEMI read in data input iii) from logical unit 8) of 


either -1/4 or -3/4 cycles is necessary. 


7.e. 


v 


sin Ae - 9.) e |b, | COs (V, 


|b, | cos Me 


ads: 1/4) if b, 2 0 
ad 3/4) if b, <i 


Similarly, an adjustment of 1/2 cycle will only be necessary for slow and 


semidiurnal main constituents if the tidal potential amplitude is negative. 


Making these changes, the combined result of a constituent cluster 


in the diurnal and terdiurnal cases is 


ja, | cos (Vi - 9;) + } As ai, COS (Vip + Way - g,) + } A 


ae Ae 
J J% Je J% J 


j Sa a ie ~ =13/4, 
where if a. <a Ww Y=) 3/4... O45 1/2 and Os 3/ 


and if as =) a 


k 


V - 1/4, ta = 0 and Ph go Sey fd. 


A further phase adjustment to satellite constituents can be made if we wish to 


ensure that their amplitudes are positive. This convention was adopted for the 


data package of Appendix 1 (variable PH read in data input iv) from logical 


unit 8). Replacing ask and as by their absolute values we now see that 


= 0 


1/2 


if both ask and a. have the same sign, 


otherwise; 


Aa 


0. 1/4 if both aay and a have the same sign, 


jer 


3/4 otherwise. 
Similarly, for the slow and semidiurnal constituents the cluster contribution 


can be rewritten as 


ee ; Viltaspads,) + ; 2 
Ja, cos (V. 93) +? Ay a5, Cosa $k 5k D5) } A la,,I cos (Vi otos 959) 


where VSR ral if ag 45 0, 
V otherwise; 
Ws = 0 ii ask and a, have the same sign, 
Lij2 otherwise; 
Aig = -1/4 Af, aso and a have the same sign 
1/4 otherwise. 


Special note should be made of the terdiurnal M3 because both it and 
its only satellite are due to third order terms in the tidal potential. Hence, 
both contribute directly through a cosine term and so behave as if they were 


second order semidiurnals. 


In order to determine the amplitude and phase of the major contributor, 
we assume that the result actually found in the analysis was f.a. cos we ~ ae 3 ui) 
where rs and u, are called the nodal modulation corrections in 
amplitude and phase respectively. To avoid a possible misunderstanding, it is 
worth mentioning here that the term nodal modulation is actually a misnomer. 
It and the symbols f and u were first used before the advent of modern computers 
to designate corrections for the moons nodal progression that were not incorpo- 
rated into the calculations of the astronomical argument for the main constituent. 
However, now the term satellite modulation is more appropriate because our 


correction is due to the presence of satellite constituents differing not only 


in the contribution of the lunar node to their astronomical argument, but also 


in the lunar and solar perigee effect. 


45 


For the purpose of calculating r and us it is assumed that the 
admittance is very nearly a constant over the frequency range within a consti- 
tuent cluster, and so Bead ie ad and tae laj,|/la;|- Lie las |/la,| 
are equal to the ratio of the tidal equilibrium amplitudes of the satellite 
to the major contributor. These ratios are latitude dependent when satellites 
of the third order are involved, necessitating the correction factors 
mentioned in 2.1.3. However, the ratios are usually small and the correction 


1S SLIgnt. 


Dropping the 'prime' notation and grouping the second and third order 
terms in one summation, the relationship between the analysed results for a 


main constituent and the actual cluster contribution is 


ead A= drs : Peek a ae | 
Fylajl cos (V;+uj-g,) = Jas] (cos (V;-9,) } a 


’ ee ee 
: K COs (V, 95405, aay) 


J 


Expanding this result and observing that it must be true for all V(t), 


the following explicit formulae are found for f and u: 


wee Mane Mek Gute 2 
feo ho ae OA } Ax rip COS (Ap + as4)) + U Ax rj, Sin (A 


2 
jk + as4)) 9 


L Ay Cak sin (Dip + 5) 
j = arctan 


ee area rpg 
Ls } Axe I 5, COS (As, ‘ “5k 


| 


For an analysis carried out over 2N + 1 consecutive observations, At 


time units apart, Ax is given by Ay = sin [(2N+]) At (ojk-95)/2] 
2N+1) sin [At 9578 /2) 


where c, is the frequency of the main contributor and 7k is that of its 
Satellite. However, ar is very nearly one even for a one year analysis and 


in the program is approximated by this value. 
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For a shallow water constituent whose frequency is calculated as 


ny C04, where o is the frequency of the yun main constituent from which 
ae Ce 
j=l 


it is derived and Cc, is the linear coefficient, the nodal modulation 


“lb 


0 
corrections for amplitude and phase are computed as f = I and 


: Noy j=} 
Uma j=l aa 


2.3.3 Final amplitude and phase results 


The result of the least squares analysis was to find for a constituent 
with frequency Oj, the optimal amplitude A, and phase o5 value for the tidal 


Signal A. cos 2m (ot ~ o;). However, due to nodal corrections, when the 
astronomical argument is calculated at the central time origin t = 0 of the 
record, we know that the actual contribution of the constituent cluster is 

¢ 


cos 2n (V. + u. - g:.). Hence, the amplitude and Greenwich phase lag of 


To ‘foal 
the constituent corresponding to frequency ue can be calculated as 


eee Nea) tan er =o G\) geieete Wee ct dice 
a ; A</ 5 and oF My us " 


2.3.4 Inferred constituents 


In accordance with previous notation, tidal signals in this section 
are assumed to be real in nature. However an alternative presentation using 
complex numbers, and the basis for the following development is given by 


Godin [1972]. 


If the length of a specific tidal record is such that certain 
important constituents will not be included directly in the analysis, provision 
is made via the data input on logical unit 4 to include these constituents 
indirectly by inferring their amplitudes and phases from neighbouring consti- 


tuents that are included. If accurate amplitude ratios and phase differences 
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are specified, inference has the effect of significantly reducing any periodic 
behaviour in the amplitudes and phases of the constituent used for the inference. 
This is due to the removal of interaction from the neighbouring inferred 
constituent. If it so happens that a constituent specified for inference is 
included directly in the analysis, the program will ignore the inference calcu- 


lations. 


The actual adjustments are as follows. Assume that the constituent 
with frequency og is to be inferred from the constituent with frequency o,, and 
that the least squares fit analysis found the latters contribution to be 
AL cos 2n(o,t - 0,°), where i, and b.. are the amplitude and phase respectively 
(o, and 0. are measured in cycles/hour and cycles respectively). Letting 

VU, be the astronomical argument + nodal modulation phase correction, 
9] be the Greenwich phase lag, 
f,; be the nodal modulation amplitude correction factor, 
and a, be the corrected amplitude, 
then from 2.3.3 we know that 
= Dies VU4. 294 
and a, = A,/f;. 
Assuming that A, and 0, are the post-inference amplitude and phase for the 
constituent with frequency oj, 
Pig = ag/az = (Ao/fo)/(A\/Ffy), 
and ¢  =/gi-g = VU; + 0) - Wp - 95 
(the latter two being data input variables R and ZETA respectively), then 


the presence of the inferred constituent in the analyzed signal yields the 


relationship: 


48 


0 0 
A, cos 2n(o,t-0; ) = Ay cos 2n(o,t-0;) + Ay we 2n(oot-0,) 


2 
A, COS 2t(o,t-0,) 1t7ry5 tiles 2n{(o5-0,)t+VU, = VU, sae 


fo 
ro A sin 2n(o,t-0,) Jn lreom 2m{(o5-0;,) t+VU, VU, ube Oo Fe 


Since the constituent with frequency o» was not chosen for inclusion 
in the least squares analysis, |o2-0,| N<RAY, where N is the record length in 
hours and RAY is the Rayleigh criterion constant (usually 1.0). Assuming in 
general that |o-o,|N is small, good approximations to 
cos 2n{(c2-0,)t + VUs-VUytc} and sin 2nf{{a5-oy)t + VU,-VU;+z} are their average 
values over the interval [-N/2, N/2], namely 
sin [mN(o2-0;)] cos [2n(VU> - VU +2)]/mN(o5-0,) and 
sin [mN(og-0,)] sin [2n(VU - VU, +z)1/mN(op-0,) respectively. 


Making these substitutions and setting 


Fo 
S = rielgelsta [mN(o-a})] sin [2n(VU, = VU, +. c)]/mN(c9-0,), 
Fo 
and Cais Pi F, sin [mN(op-0))] cos[2n(VU, - VU, +¢)]/mN(ap-0)), 


we obtain 
: 0 
a cos 2n(o,t - 0; ) = C cos 2n(oyt - 0) - S sin 2n(o,t - Oy). 
Ay 


£xpanding and regrouping this result yields 


0 
A 


1 | 0 
cos aot ge om 2m, -C cos 2n0,- S sin Agri 
| 


0 


Ay ( 
sin anit e Sin 2nfy +C sin 2n0, - S cos at). 
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Now since this relationship must hold for all t, both terms in 


brackets are equal to zero. Hence 


0 
0 
a cos 2m, = C cos 2nh, + S sin 2n0,, 
Aj 
0 
0 
> Sin 219; = C sin 2b, - S cos 2nO,, 
Ay 


and so 
Ve Cay Marware yg 
0; = 1 + [arctan (S/C)]/2n. 
The relative phase and amplitude of the inferred constituent are then calculated 
as 
d2 = VUy - Wy + 0) -c, 
and As = ry A, (fo/f,). 


ot 


Soc 


aS) 
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USE OF THE TIDAL HEIGHTS PREDICTION COMPUTER PROGRAM 
General Description 

This program produces tidal height values at a given location 
for a specified period of time. Amplitudes and Greenwich phase lags 
of the tidal constituents to be used in the prediction are required 
as input and either equally spaced heights or all the high and low 


values can be produced. 


Routines Required 


LOU MAIN GD eee ee eee reads in tidal station and time period 
information, amplitudes and Greenwich phases 
of constituents to be used in the prediction, 
and calculates the desired tidal heights. 


2): ASTRO Gate etre tereets ee reads the standard constituent data package 
and calculates the frequencies, astronomical 
arguments, and nodal corrections for al] 


constituents. 
3) TPUTT eaiorse, «arene saan beenvciere controls the output for high-low predictions. 
Ai) THRU aa ecc cee ates een one controls the output for equally spaced 
predictions. 
Bly CONV era < cated ree ateeetene returns the consecutive day number from a 


specific origin for any given date in the 
range 1901 to 1999, and vice versa. 


Data Input 

All input data required by the tidal heights prediction program is 
from logical unit 8. A sample set is given in Appendix 4. Although 
data types i), ii) and iii) are identical to types ii), iii) and iv) 
expected on logical unit 8 by the analysis program, for completeness they 


are repeated here. 


i) 


5] 


Two cards specifying values for the astronomical arguments SO, HO, PO 


>’ 


ENPO, PPO, DS, DH, DP, DNP, DPP in the format (5F13.10). 


SOQ. = mean longitude of the moon (cycles) at the reference time 
origin. 

HO = mean longitude of the sun (cycles) at the reference time 
origin. 

PO = mean longitude of the lunar perigee (cycles) at the reference 


time origin. 


ENPO = negative of the mean longitude of the ascending node (cycles) 
at the reference time origin. 


PPO = mean longitude of the solar perigee (perihelion) at the 
reference time origin. 


DS, DH, DP, DNP, DPP are their respective rates of change over a 


365 day period at the reference time origin. 


At least one card for all the main tidal constituents specifying 
their Doodson numbers and phase shift along with as many cards 

as are necessary for the satellite Eee init s:, The first card 
for each such constituent is in the format (6X,A5,1X,613,F5.2,14) 


and contains the following information: 


KON constituent name, 
II],JJ,KK,LL,MM,NN = the six Doodson numbers for KON, 


SEMI 


the phase correction for KON, 


NG the number of satellite constituents. 


A blank card terminates this data type. 
If NU > 0, information on the satellite constituents follows, three 
satellites per card, in the format (11X,3(313,F4.2,F7.4,1X,17,1X)). 


For each satellite the values read are: 
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LDEL, MDEL, NDEL = the last three Doodson numbers of the main 
constituent subtracted from the last three Doodson 
numbers of the satellite constituent; 


PH = the phase correction of the satellite constituent 
relative to the phase of the main constituent; 


EE = the amplitude ratio of the satellite tidal 
potential to that of the main constituent; 


IR = ] if the amplitude ratio has to be multiplied 
by the latitude correction factor for diurnal 
constituents, 


= 2 if the amplitude ratio has to be multiplied 
by the latitude correction factor for semi- 
diurnal constituents, 


otherwise if no correction is required to the 

amplitude ratio. 

iii) One card specifying each of the shallow water constituents and the 
main constituents from which they are derived. The format is 


(6X ,A5,11,2X,4(F5.2,A5,5X)) and the respective values read are: 


KON = the name of the shallow water constituent, 


NJ = the number of main constituents from which it 
is derived, 


COEF, KONCO = combination number and name of these main 
constituents. 


The end of these shallow water constituents is denoted by a blank card. 


iv) One card with the tidal station information ISTN, (NA(J),J=1,4), 
ITZONE, LAD, LAM, LOD, LOM in the format: 
(5X5 1491X, 3A6 ,A4SASAIK, LARRX LL 2M LS ENS LZ) 3 
ISTN = the station number. 
(NA(J),J=1,4) 


the station name. 


ITZONE = the time zone reference for the 'Greenwich' phases. 


LAD, LAM the station latitude in degrees and minutes. 


LOD, LOM 


the station longitude in degrees and minutes. 


v) 


vi) 


03 


One card for each constituent to be included in the prediction, 

with the constituent name (KON), amplitude (AMP) and phase lag 

(G) in the format (5X,A5,28X,F8.4,F7.2). (This format is compatible 
with the analysis program results produced on output device 2). The 
phase lag units should be degrees (measured in time zone ITZONE) 
while the units of the predicted tidal heights will be the same as 


those of the input amplitudes. The last constituent is followed 


by a blank card. 


One card containing the following information on the period and 


type of prediction desired. The format is (313,1X,313,1X,A4,F9.5). 


IDYO, IMOO, IYRO = the first day, month, and year of the 
prediction period. 


IDYE, IMOE, IYRE = the last day, month and year of the 
prediction period. 


TTYPE = 'EQUI' if equally spaced predictions are desired, 


= "EXTR' if all the high and low tide times and 
heights are desired.» 


DT = the time spacing of the predicted values if 
Leer "FOUL 


= the time step increment used to initially 
bracket a high or low value if ITYPE = 'EXTR'. 

Equally spaced predictions begin at DT hours on the first day and 

extend to 2400 hours (assuming 24 is a multiple of DT) of the 

last day. When ITYPE = 'EXTR', Godin and Taylor [1973] recommend 

using the following values for DT: 3 hours for a semi-diurnal tide, 


6 hours for a diurnal tide, and 0.5 hours for a mixed tide. 


Type vi) data may be repeated any number of times. One blank card 
following a type vi) record will return the program to type iv) 


input, while two blank cards will end the program execution. 


3.4 
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Output 


Two logical units are used for the output of results in the tidal 
heights prediction program. Device number 6 is the line printer and 10 
is a data file. Both equally spaced and high-low predictions are put 
onto both devices with the same format. However the line printer also 
records the station name and location along with the amplitudes and 
phase lags of the constituents used in the prediction. Appendix 5 


lists device 10 output resulting from the input of Appendix 4. 


When daily high-low values are desired, the date, station number 
and a series of up to six heights and occurrence times are listed 
per record. Each record begins with the variable HL whose value is 0 
if the first height for that day is a high (i.e. larger than the second 
height) and 1 if the first height is a low. If there are less than six 
high-low values for a day, they are padded up to six with the values 
9999 and 99.9 for the times and heights respectively. On device 10, 
the format used for the variables HL, the station number, the day, 
month, year, and six times and heights is 

(1X, Ul Seales D2eG i bebe i) 

When equally spaced heights are requested, 8 values are listed on 
each record preceded by the station number, the time, day, month and 
year of the first value, and followed by the time increment between 
heights. On device number 10, the format for these variables is; 


GUAM T Ae as. 21 aror6.2 2A) 


) 


a0 


Program Conversion, Modifications, Storage and Dimension Guidelines 


The source program and constituent data package described in 


this manual were tested on the UNIVAC 1106 at Patricia Bay Institute 


of Ocean Sciences. Although the program was written in basic ASCII 


FORTRAN, there may be some changes needed before the program and data 


package can be used on other installations. These may include: 


i) 


ii) 


iii) 


converting the character code from EBCDIC to BCD. 


deleting some or all the REAL*8 declaration statements. Presently 
these are used only to permit alphanumeric strings longer than 4 
characters (which is the ASCII FORTRAN single precision word 
length) to be read (eg. constituent names are read in A5 format). 
CDC installations should not require these double precision 
variables because there are no alphanumeric strings longer than 


6 characters (which is the CDC word length) in the program. 


Switching the Chebyshev iteration formula to double precision 
in order to maintain accuracy. Tests performed on the UNIVAC 
demonstrate that this change is not necessary, however on less 
accurate installations, it may be (UNIVAC real variable single 
precision has 8.1 significant digits). As a test for machine 
accuracy, I suggest using the input of Appendix 4 and comparing 
the results with those listed in Appendix 5. In the event that 
more precision is needed, the following variables in the main 
program should be declared REAL*8: CH,CHA, CHB, CHM, CHP, SUM 
TWOG, BIWOGey HGR 2, <ZAy- ZBy°ZP 502M, -DTPY; 

and all related assignments and library function calls should be 
altered accordingly (e.g. '=0.0' should become '0.D0' and SQRT 
should be DSQRT). 
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iv. altering the variable list structure for ENTRY statements and 


references to them. 


v) changing some or all of the logical unit (file reference) numbers 


from their present values in order to conform with local machine 


rescricLionss 


The program in its present form requires approximately 3000 and 


15000 UNIVAC words for the storage of its instructions and arrays 


respectively. 


As with the analysis program, changing the number or type 


of constituents in the standard data package may require alteration to 


the dimensions of some arrays. Restrictions on the minimum dimension of 


all arrays are now given. 


Let 
MTAB 
M 
MCON 
MSAT 
MSHAL 
and NITER 


be the total number of possible constituents contained in 
the data package (presently 146), 


be the number of constituents to be included in the 
prediction, 


be the number of main constituents in the standard data 
package (presently 45), 


be the sum of the number of satellites for these main 
constituents and the number of main constituents with no 
satellites (presently 162 + 8) 


be the sum for all shallow water constituents of the number 
of main constituents from which each is derived (presently 
25 lai 


be the number of iterations required to reduce the time 
interval within which it is known that a high or low tide 
exists, to a desired length (with the largest initial 
interval size of 6 hours and a 6 minute final interval, 
NEVER: 1s 16). 


7, 


Then in the main program, arrays KONTAB, SIGTAB, V, U, and F should 
have minimum dimension MTAB; KON, SIG, AMP, G, INDX, TWOC, CH, CHP, CHA, 
CHB, CHM, ANGO and AMPNC should have minimum dimension M; and the two 
dimensional array BTWDC should have a minimum dimension of M by NITER. 
Array COSINE which stores pre-calculated cosine function values over the 
range of 0° to 360° and is used as a look-up table, presently has 2002 
elements. 

In subroutine ASTRO, the arrays KON, FREQ, V, U, F and NJ should 
have minimum dimension MTAB; arrays II, JJ, KK, LL, MM, NN and SEMI 
should have minimum dimension MCON; arrays EE, LDEL, MDEL, NDEL, IR and 
PH should have minimum dimension MSAT; and arrays KONCO and COEFF should 


have minimum dimension MSHAL. 


In subroutine PUT, the dimensions of arrays HGTK and ITIME should 
be at least as large as the maximum number of high and low values per 


day (this is presently assumed to be 9). 


In subroutine HPUT, the dimension of array H should be at least 
equal to the number of equally spaced tidal height values per output 


record of logical unit 10 or 6 (presently, this is 8). 


In subroutine CDAY, both arrays NDM and NDP should have dimension 12. 
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4. Tidal Heights Prediction Program Details 


4.1. Problem Formulation and the Equally Spaced Predictions Method 
The tidal height, h(t), at a particular station may be represented 


by the harmonic summation (see section 2.3.3) 


h(t) =") F.CtYA, cos(2n(M(t) + welt) eige)) os 1 
(t) = JE) Fy (t)A; cos(2n(v;(t) + uj(t) ~ 95) (1) 
where Aas a - the amplitude and phase lag of constituent j, 
f(t), u(t) = the nodal modulation amplitude and phase correction 
factors for constituent j, 
and V(t) = the astronomical argument for constituent j. 


Expanding V(t) as in section 2.3.1 and using the first order Taylor 
approximations for the astronomical arguments as in section 2.1.1, V(t) can 


be re-expressed as 


V(t) = ix(t) pegS(e)meekH(t) eee (t) SpemNa (tyr onk st 


it(to) + JS(to) + kH(to) + 2P(to) + mN'(t.) + nP' (to) 


+ (t-te Lic(t)+jS(t)+kH(t)+2P(t)+mN' (t)+nP'(t) Jy 


Ste) tb aGtee Lovo. 


where t. is the reference time origin of 000 ET January 1 1976, and o is the 
constituent frequency at this time origin. It follows from this result that 
V(t>) = V(t,) + (t. - t,)o for arbitrary times t,, t,, and so V(t) can be 


replaced in (1) by Ney) ls ae tio, for some convenient time t,. 
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From section 2.3.2 it is seen that f(t) and u(t) are time dependent 
only through the As, (t) variable. Since satellites differ from main 


constituents in only the last three Doodson numbers (see section 2.1.3), 


ul 
=< 
~~ 
ct 
~~ 
\ 
=< 
| 
ct 
~~ 


A(t) jk 


ARP(t) + AmN'(t) + AnP'(t) . 


Using the first order Taylor approximations for P,N', and P', it follows that 
over a time period [t,,t.] the change in Aa, (t) is 


ds, (te) - As, (ti) = Ag(P(t.)-P(t,))+am(N'(t,)-N'(t,))tan(P'(t,)-P'(t,)) 


= (to-ty) GefaeP(t) + amN'(t) + Ae alene 


= (to-ty )(o5,-05). 


Since [P(t)+N'(t)+P'(t)] 4 is .16668884 cycles/365 days and |Ag|,|Am],|An| 


are always less than or equal to 4, if |t,-t,|<!6 days, |Asy (ta)-a5, (ti) [5 .03 cycles. 
This small variation in As, (t) leads to a similar behaviour in cos(A,,(t)) and 
sin(4,,(t)), and hence f(t) and u(t). Thus only a small loss in accuracy but a 
considerable calculation time saving will result if f(t) and u(t) are 
approximated by a constant value throughout the period of a month. Consequently 
f(t) and u(t) are assumed to equal their value at 000 hr of the 16th day of the 
month for the entire monthly period; and for convenience, V(t) is set to 


V(tig) + (t-ty¢)o, where ty), is this same time. 
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The procedure for calculating a series of tidal heights is then as 
follows. Since the tidal prediction data package does not contain constituent 
frequencies, they must be calculated via the astronomical variable derivatives 
and the constituent Doodson numbers. f, u and V values are then calculated 
for the 16th day of the first month of the desired prediction period and, 
as required, for subsequent months. Tidal heights for the desired values of 


t can then be calculated as 


m 
h(t) = bi F (tig )Ajcos(2m(V 5(ty6)#(t-tyg)o stu; (ti6)-95)). (2) 


In order to avoid calling a trigonometric library function for each new value 
of t, when a sequence of equally spaced heights are required, the following 


Chebyshev iteration formula is used for each constituent contribution, 
Fintiver2-cos(cAt tine -sbin-1) . (3) 


where f(n) = cos(noAt) or sin(nost). 
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4.2 The High and Low Tide Prediction Method 


The material presented here is taken from Godin and Taylor [1973]. 
In 4.1 we saw that the tidal height at a given location can be 
represented by the harmonic sum 


h(t) = D521 fj(to)Ajcos(2m(V;(to) + (t-to)o, + u(te) - g,)) (1) 


where: 


Sete oe = the amplitude, phase lag and frequency of constituent j, 
f.(t.),u.(t.) = the nodal modulation amplitude and phase correction 


factors for constituent j at the time origin to, 


<= 
Ce 
ct 
° 
~ 
! 


= the astronomical argument for constituent j at the time 
OF OI bot: 
Letting D(t) be the derivative of h(t), i.e. 


D(t) ee ey F5(to) A j@na;sin(2n(V (to) a (t-to)o. mE u(t.) 95)) 9 


the high-low tide prediction method uses the following calculus results. If 


D(t) is a continuous function on the interval [t,,t»] and ty is a point in this 


interval, then: 
i D(t,) = 0 if and only if ty is an extreme point, or saddle point,! or 
h(t) is constant in the neighbourhood of ty 
ii) if D(t,) and D(t,) have opposite signs then there exists a ty in 


Lr casliod) with D(t,) 


1 An example of a saddle point is x = 0 for the function f(x) = x9. 
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Now for computational purposes we can assume that saddle points do not 
exist. That is to say, due to accuracy limitations of the computer, a zero 
derivative will be approximated by a number with a very small absolute value, 
and thus perturb a saddle point so that it becomes either a maximum and 
minimum, or a near saddle point (in the neighbourhood of a ‘near saddle point’, 
the derivative is of constant sign and almost assumes the value zero). And 
since, from its definition, we can reasonably assume that h(t) is not constant 
over any arbitrarily small interval, the continuity of D(t) everywhere implies 
that an interval [t,,t.] with D(t,) and D(t,) reaian opposite signs contains an 


extremum. 


However this result alone is not sufficient to guarantee the location 
of all extrema because, it does not eliminate the possibility of having more 
than one extremum in an interval whose endpoints have different signs, nor does 
it imply that if the endpoints have the same derivative sign there is no 
extremum in the interval. In order to ensure these conditions and thus be 
assured of bracketing all extreme values, it is necessary that a minumum interval 
size be specified in which we can assume that there exists at most one high or 


low tide. 


Clearly the interval size At will be dependent upon the nature of the 
tide at a particular station. The time between successive high and low waters 
for predominantly semi-diurnal and diurnal tides is approximately 6 and 12 hours 
respectively. However if the tide is mixed, the extremal pattern is more 
complicated. FIGURE 2 shows the water level at Victoria, British Columbia 
between July 24 and 31, 1976. It is a mixed tide where the shorter period 
fluctuations override the major diurnal oscillations with a continuous shift in 


their position and amplitude. 
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One characterization of the tide may be obtained by calculating the 
ratio of the amplitudes of the major harmonic constituents M2, S2, Kl and Ol. 


This value is called the form number (Dietrich [1963]) and is defined precisely 


as 
oy SKA aD 
Meee 
The tide is then said to be 
i) semi-diurnal if Ose Fee coy 
ii) mixed if 125) <thee so O08 
iii) diurnal Th Pies s00r- 


For Victorias Preval. 
In accordance with this determination Godin suggests the following 
maximum time interval values in which it can be assumed that there exists at 


most one extremum: 


i) At = 3 hours for semi-diurnal tide, 


ii) At = .5 hours for mixed tide, 


iti At, = "6 hours tor ciurnal tide. 


Although in fact, a mixed tide may have extrema closer than .5 hours, he feels 


that for practical purposes it is sufficient to note just one of them. 


With these values of At we can then bracket all extrema by moving 
forward in time with steps of size At, and comparing signs of the interval 
endpoints. Once such upper and lower bounds have been found, the extreme 
point can be located exactly by any one of a number of search techniques. 


Because it requires a minimal amount of time, the one chosen is Bolzano's 
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method of bisection coupled with linear interpolation. Although the bisection 
method does not take the minimal number of iterations when compared to more 
sophisticated search techniques, it is able to make significant time savings 
by computing new sine function values as a linear combination of old ones, and 


thus, unlike the other methods, avoid calls to the fortran library function SIN. 


In more detail, the search algorithm for an extremum is then as follows: 
i) move forward in time from the origin, or the last extremum, in steps 

of At until either a change in sign exists between the derivative 

values at the end points of the interval (t,.t,); or t, extends 

beyond the desired prediction period. Each constituent contribution 

in the summation D(t) is evaluated by the Chebyshev iteration 

formula (3) of 4.1. When an interval containing an extremum is 


located, set k = 1 and proceed to ii). 


] 
ii) calculate tL = t, + ok At, and for each constituent in the sum, 


evaluate D(t, ) by using the formula 


sin(t, ) = (sin(t,) + sin(t,))/(2 cos (ey At)). 


If |D(t,) | SH OFS puseteDity) sehOn? 


Q) 


iii) re-assign whichever of te or th has the same derivative sign as D(t,), 


by t If the new interval length tt. is less than .1 hours, 


Kk 


proceed to iv). Otherwise set k = k+] and return to ji). 


iv) use the following linear interpolation formula to find the extremum 


fer oe t D(t,)(t,-t.)/(D(t 


: ech 


a 
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and evaluate h(t.) via (1). For each constituent term in this sum, 
obtain the function value by using a pre-calculated stored table of 
2002 cosine values with arguments in the range of 0° to 360°. Return 


tOeP TY): 


FIGURE 3 illustrates an example of the sequence of steps involved in 
the search for an extreme value. It is easily calculated that the number of 
iterations required to reduce the bracketing interval from At to .1 hour is 


6 for diurnal tides, 3 for mixed tides, and 5 for semi-diurnal tides. 


Arrows in FIGURE 2 indicate the extrema predicted for Victoria using 
the technique just described; the shaft of the arrow locates the time abscissa 
while the tip ends at the predicted height. The predicted hourly heights and 


the times and heights of all extrema are listed in Appendix 5. 
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An example of the sequence of steps involved in locating a zero 
ty of the derivative D(t). The sign of D(t) at the various points tested 


is denoted by + or -. After a step At, the sign has changed; by a 
retrogression of At, the sign has reverted to plus, forcing a forward step 
of 4At where the sign is still unchanged. Two further forward steps of 
Ikat and 4AgAt locate the minimum width interval (tot, ) over which the 
position of ty is determined by linear interpolation from the values of 


OCLs at t. and tad: 
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5. Consistency of the Analysis and Prediction Programs 


Although consistency between the tidal heights analysis program and 
the tidal heights prediction program was a major objective in their revision, 
they do have one difference. In particular, if a pseudo-tidal record were 
generated by the prediction program and analysed using the same constituents, 
the amplitude and phase results given by the analysis program would not be 


identical to those used as input for the prediction program. 


In a small part, this discrepancy is due to round-off accumulated 
during the calculations. However a test performed on the UNIVAC 1106 at 
Patricia Bay with a two month period of synthesized hourly heights indicates 
that such errors occur no sooner than the fourth digit. The remainder of the 
difference (which is, at worst, in the third digit) can be attributed to 
different approximating assumptions for the calculation of f and u, the nodal 
modulation amplitude and phase correction factors. Whereas the prediction 
program calculates these values at the 16th day of each month in the desired 
time period and keeps them constant aycauretia: the entire month, the analyses 
program assumes them to be constant over the entire analyses period and equal 


to their true values at the central hour of that period. 


It is important to note, though, that significantly different results 
can be expected in a similar test run if there is at least one more constituent 
used in the synthesis than analysis. This is because the least squares fit 
technique will adjust the amplitudes and phases of constituents included in 


the analysis to partially account for contributions due to constituents 
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included in the synthesis but not the analysis. In fact this will occur even 
if the extra constituents are inferred (e.g. Pl is included in the synthesis 
and in the analysis via inference from Kl) because of smal] inaccuracies in 
the approximating inference assumptions. However, except for round-off errors 
and the slightly different f and u values, having more constituents in the 


analysis than the synthesis will not affect the results. 
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